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Modeling  for  Precise  Terrestrial  Inertial  Navigation".  This  report 
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to  the  School  of  Engineering,  Air  Force  Institute  of  Technology, 
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I.  Introduction 


Inertial  navigation  has  evolved  to  the  point  that 
the  traditional  gravity  model  is  a  principle  error  source 
in  advanced,  precise  systems.  Future  autonomous  strategic 
systems,  such  as  the  cruise  missile,  will  likely  require 
a  more  accurate  gravity  model.  Mission  success  depends 
on  navigation  system  accuracy  not  gravity  model  accuracy 
per  se.  So  an  improved  gravity  model  should  be  judged  by 
system-level  performance  expressed,  for  example,  as  circular- 
error-  probable  (CEP).  The  objective  of  this  research  is 
to  develop  a  new  computational  technique  whereby  alterna¬ 
tive  gravity  models  can  be  compared  by  expected  navigation 
system  accuracy  on  realistic  mission  scenarios. 

This  section  presents  a  discussion  of  the  background 
and  motivation  for  the  research.  The  role  of  the  gravity 
model  and  the  navigation  system  errors  induced  by  this 
model  are  interpreted  for  this  study.  The  need  for  a 
statistical  evaluation  is  explained,  and  the  present  methods 
for  performing  such  analyses  are  discussed.  The  short¬ 
comings  of  these  present  statistical  evaluation  methods 
are  given  as  motivation  for  the  research. 

1.1  The  Role  of  the  Gravity  Model  in  Inertial  Navigation 

The  gravity  model  in  an  inertial  navigation  system 
is  a  subsystem  or  component  providing  one  part  of  the 
information  required  for  navigation  position,  velocity  and 
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attitude  estimates  (see  Figure  1).  The  exact  role  of 
the  gravity  model,  the  nature  of  modeling  errors,  and  the 
effect  of  such  errors  on  navigation  estimates  must  be  under¬ 
stood  to  appreciate  the  impetus  for  gravity  model  improve¬ 
ments  . 

An  inertial  navigation  system  senses  the  vehicle 
dynamics  and  applies  the  laws  of  physics  to  estimate  posi¬ 
tion,  velocity  and  attitude.  The  inertial  navigation 
system  performs  two  fundamental  tasks.  First,  the  vehicle 
rotational  velocity  ^  is  measured  and  used  to  track  attitude. 
Second,  position  r  and  velocity  v  are  computed  by  inte¬ 
grating  an  estimate  of  acceleration  a. 

The  total  vehicle  acceleration  cannot  be  sensed, 
since  the  gravity  field  acts  on  both  the  vehicle  and  the 
accelerometers.  The  total  acceleration  is  given  by 

a  =  f  +  G  (1) 

where  f  is  the  specific  force  (force  per  unit  mass)  from 
contact  forces  such  as  aerodynamic  drag  and  G  is  gravita¬ 
tional  acceleration  due  to  mass  attraction  acting  on  the 
vehicle.  The  estimate  of  acceleration  is  made  by 

a  =  f  +  Gm(r)  (la) 

where  "A"  indicates  navigation  estimates,  indicates 

measured  data,  and  G  (•)  is  the  gravitation  model. 

The  gravity  model  term  in  (la)  introduces  error  by 

a.  Modeling  errors,  and 

b.  Evaluation  errors 


The  evaluation  errors  come  from  evaluating  gravitation  at 
an  erroneous  position  r.  This  evaluation  error  is  respon¬ 
sible  for  the  principle  navigation  system  error  fnode  -  the 
Schuler  mode.  The  evaluation  error  is  not  an  error  induced 
by  the  model  and  is  therefore  considered  part  of  the  in- 
ertial  navigation  error  propagation.  The  modeling  error 
can  be  st^ifed  as  a  gravity  disturbance 

6g(r)  =  Gm(r)  -  G(r)  (2) 

given  pei*fect  position  information.  This  error  of  (2) 
introduced  by  (la)  creates  the  navigation  system  errors 
which  are  the  subject  of  this  work. 

1.2  The  Traditional  Gravity  Model 

The  gravity  model  used  in  most  present  operational 
inertial  navigators  is  based  on  a  gravity  field  perpendic¬ 
ular  ("normal"  field)  to  an  ellipsoid  which  approximates 
the  mean  sea-level  equipotential  surface  called  the  geoid. 

In  some  instances,  zonal  spherical  harmonic  or  spheroidal 
equipotential  surface  models  are  used  as  a  basis  for  the 
gravity  model.  The  associated  simple  gravity  models  are 
within  400  mgal*  of  measured  gravity  values. 

The  use  of  such  simple  models  has  been  prevalent  for 
two  reasons.  First,  computational  resources  are  limited 
for  the  on-line  application.  Second,  the  presence  of  other 
system  error  sources  has,  in  the  past,  reduced  the  incentive 


*  1  mgal  =  10“3  galileo  =  10~3cm/sec2  =  10~3m/sec2  2  lpg 
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for  a  more  complex  and  accurate  model.  Kay ton  [Ref  l]  in 
i960  suggested  that  these  simple  models  are  sufficient 
for  system  with  accelerometer  uncertainty  greater  than 
20  mgal  —  an  estimate  at  that  time  of  gravity  disturbance 
standard  deviation. 


1.3  Impetus  for  Model  Improvement 

Since  the  late  1940' s,  inertial  navigation  design  has 
progressed  to  a  fine  art.  Refinements  have  been  concen¬ 
trated  principally  in  inertial  instrument  design  and  cali¬ 
bration  with  other  basic  concepts  essentially  intact.  The 
error  levels  of  periodically  recalibrated  inertial  components 
have  reached  the  point  that  the  traditional  gravity  model 
is  now  a  principle  error  source  in  overall  system  accuracy. 

It  is  easy  to  understand  why  a  detailed  new  model  might  be 
used  in  a  flight  test  environment  where  data  purity  is 
important.  However,  one  could  question  the  need  on-line 
for  model  refinements  since  the  errors  induced  are  not 
unacceptable  for  most  navigation  applications.  The  impetus 
comes  from  the  potential  military  application  —  where  the 
impetus  for  inertial  navigation  originated.  In  the  autono¬ 
mous  delivery  of  weapons,  any  error  significantly  diminishes 
weapon  systems 's  effectiveness,  so  these  gravity  induced 
errors  cannot  be  ignored.  The  self-contained  nature  of 
inertial  navigation  virtually  assures  a  continued  military 
application  even  with  advanced  radiometric  navigation  sys¬ 
tems  available.  The  evolution  of  strategic  forces  to 


include  cruise  missile  concepts  provides  the  prime  motiva¬ 
tion  for  the  increased  emphasis  on  refining  the  traditional 
gravity  modeling  used  in  on-line  inertial  navigation. 

1.4  The  Need  for  Model  Performance  Evaluation 

The  need  for  an  improved  gravity  model  may  be  met 
using  existing  data  base  and  on-line  modeling  methods  CRef 
2].  While  future  research  may  refine  and  make  them  more 
efficient,  several  methods  of  computing  gravity  exist  today. 
These  methods  are  primarily  based  on  Green’s  and  Stoke’ s 
theorems  [Ref  3]  which  require  data  over  a  closed  surface 
(approximately)  encompassing  all  earth  mass.  Heretofore, 
data  limitations  over  vast  ocean  expanses  have  diminished 
the  accuracy  of  these  models.  Recent  satellite  altimetry 
measurements  give  an  estimate  of  the  anomalous  potential* 
over  much  of  these  regions.  These  new  data  can  be  combined 
with  existing  gravity  data  using  heterogeneous  data  proces¬ 
sing  techniques  [Ref  43.  This  increased  data  base  can  then 
be  used  to  identify  parameters  of  a  more  precise  model  for 
on-line  application. 

The  gravity  models  which  might  be  used  on-line  have 
many  forms.  The  canonical  form  is  the  Legendre  function 
expansion  based  on  a  Fourier  series  representation  of  the 
geoid.  This  spherical  harmonic  series  can  theoretically  be 
expanded  to  any  degree  and  order  required  to  yield  a  residual 


*  See  Reference  3  and  Appendix  A  for  definitions  of  anomal¬ 
ous  field  quantities. 
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field  below  some  arbitrary  magnitude  33  •  Measurement 

errors  and  computational  costs  argue  against  carrying  such 
modeling  to  the  extreme,  however.  This  form  of  globally 
applicable  functional  approximation  is  countered  by  other 
methods  which  seek  to  describe  gravity  in  a  smaller  region 
using  interpolating  functions  [Ref  53-  The  rationale  is 
that  fewer  coeff icients ,  hence  fewer  calculations,  may  be 
required  for  the  same  accuracy  within  a  restricted  region. 
The  point  mass  grid  of  the  MINUTEMAN  launch  region  gravity 
model  [Ref  6]  is  a  compromise  between  these  extremes. 

While  this  model  is  a  globally  valid  approximation,  it  has 
more  detail  and  greater  accuracy  in  the  critical  launch 
region.  The  Poisson  integral,  Stoke 's  integral,  and  the 
coating  method  [Ref  2:39-5^3  give  a  more  direct  link  between 
measured  data  on  one  hand  and  gravity  estimates  on  the 
other.  Approximations  of  field  theory  integrals  consti¬ 
tute  alternate  forms  of  gravity  modeling. 

1.5  Present  Methods  for  Evaluating  Alternative  Gravity 

Models _ 

An  abundance  of  modeling  methods  exists  then.  What 
is  lacking  is  a  consistent  way  of  comparing  models  and 
specifically  evaluating  the  error  contribution  to  inertial 
navigation  performance  corresponding  to  each  candidate 
model.  There  are  two  fundamentally  different  methods  of 
calculating  a  system  accuracy  measure  1  deterministic  and 
statistical.  In  a  deterministic  analysis  a  complex  gravity 


model  is  used  as  the  "truth"  model.  The  design  mission 
is  simulated,  and  the  resulting  inertial  navigation  esti¬ 
mates  are  compared  to  the  true  position,  velocity  and  atti¬ 
tude.  This  deterministic  study  yields  error  as  a  function 
of  time.  The  resulting  error  profile  is  only  valid  for 
that  particular  mission  and,  hence,  should  be  considered 
one  (simulated)  sample  out  of  the  set  of  possible  missions. 
Since  the  truth  model  is  limited  in  scope  (it  still  has 
errors  of  omission)  and  has  an  inherently  inaccurate  data 
base  (errors  of  commission  cannot  be  avoided) ,  some  resid¬ 
ual  uncertainty  remains  in  even  this  result. 

The  statistical  approach,  on  the  other  hand,  frankly 
admits  these  model  uncertainties  at  the  outset  and  proposes 
to  project  these  model  errors  in  terms  of  navigation  per¬ 
formance.  To  this  end,  the  statistics  of  the  residual 
field  (correlation  function)  coupled  with  a  model  of  either 
the  inertial  navigator  per  se  or  of  the  error  propagation 
in  the  inertial  navigator  are  used  to  produce  estimates  of 
inertial  navigation  error  statistics.  The  general  strategy 
of  statistical  analysis  is  presented  in  Figure  2. 

1.6  Trajectory  Models 

The  trajectory  m'odel  includes  all  environmental, 
dynamic  and  ambient  condition  data  required  for  analyses. 
Usually,  position,  velocity,  specific  force,  and  attitude 
are  sufficient.  Other  conditions  such  as  temperature 
may  influence  the  operation  of  the  inertial  system  enough 
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to  warrant  modeling.  For  the  basic  problem  of  providing 
trajectory  information,  two  clear  methods  exist. 


For  relatively  simple  trajectories,  closed-form 
mathematical  expressions  can  be  derived  and  offer  a  compu¬ 
tationally  efficient  manner  of  providing  the  required 
data.  Examples  of  this  type  of  model  are  used  in  Section 
IV  studies. 

While  complex  trajectories  could  theoretically  be 
decomposed  into  a  finite  number  of  simple  segments,  deriving 
closed-form  expressions  becomes  an  onerous  burden.  Computer 
approximation  based  on  the  kinematic  relationships  of  such 
segments  can  be  programmed  and  offers  a  flexible  method  of 
generating  the  design  mission  data.  The  PROFGEN  [Ref  7], 
profile  generator,  program  developed  by  Musick  is  an  excel¬ 
lent  example  of  this  approach.  The  program  accepts  mission 
phase  specifications  as  inputs.  The  results  for  contiguous 
segments  are  abutted,  and  the  overall  mission  trajectory 
is  the  output. 


Mission 

Phase 

Specifications 


b 


PROFGEN 


r(t) 
v(t) 
etc . 


Adequate  methods  exist  to  furnish  the  trajectory 
information  required  in  the  statistical  analysis.  The  great 
flexibility  of  these  methods  should  handle  any  conceivable 
navigation  mission.  Examples  are  presented  in  later  sec¬ 
tions  using  both  closed-form  trajectories  and  PROFGEN 
approximations . 


The  inertial  navigation  system  has  been  modeled  by 
two  methods  which  should  be  adequate  for  any  statistical 
analysis;  whole-valued  simulation  and  error  propagation 
model.  The  navigation  system  model  used  may  be  dictated 
by  the  analysis  method,  but  analyses  have  been  conducted 
with  both  of  these  model  types. 

The  whole-valued  simulation  is  a  direct  simulation 
of  the  inertial  navigation  algorithm.  Errors  are  simu¬ 
lated  in  the  input  data  to  the  algorithm  and  are  observed 
in  the  output  navigation  estimates.  In  modern  navigation 
systems,  the  algorithm  is  already  in  the  form  of  a  digital 
computer  program,  so  the  operational  algorithm  may  be  used 
directly.  More  generally,  the  double  integration  of 
acceleration  may  be  used  as  a  generic  navigation  system 
model  with  no  direct  link  to  an  existing  system. 

The  second  method  is  simply  a  first-order  approxima¬ 
tion  of  the  propagation  of  errors  through  the  inertial 
navigation  system.  Only  the  errors  and  not  the  whole¬ 
valued  estimates  are  modeled.  Britting  has  shown  that 
this  first-order  error  model  can  be  formulated  in  a  general 
manner  CRef  8j.  This  formulation  can  be  cast  in  the  form 
of  a  first-order,  vector  differential  equation  CRef  9» 
Equation  3-1,  p22] 


where 


x(t)  is  an  n-.state  mathematical  vector  of  navigational 
errors  in  position,  velocity,  and  attitude  esti¬ 
mates 

F(t)  is  an  nXn  error  propagation  matrix  relating  error 
rate  to  present  state; 

u(t)  is  an  m- dimensional  vector  of  gravity  disturbance 
quantities;  and 

G(t)  is  an  nXm  distribution  matrix. 

Note  that  for  this  work,  u(t)  contains  only  those  terms 
associated  with  geodetic  errors.  For  example,  u(t)  =  6g(t). 
Other  system  error  sources  such  as  accelerometer  bias  and 
gyroscopic  drift  are  assumed  to  be  zero,  so  the  results 
obtained  here  will  reflect  the  effects  of  geodetic  errors 
alone. 

The  whole-valued  simulation  contains  nonlinearities 
which  the  error  propagation  model  loses  by  the  first-order 
expansion.  The  error  propagation  model  has,  however,  been 
widely  accepted  since  it  accurately  models  navigation  sys¬ 
tem  behavior  observed  in  flight  tests.  Either  of  these 
methods,  then,  will  sufficiently  model  the  inertial  naviga¬ 
tion  errors  to  support  the  statistical  analysis. 

1.8  Gravity  Disturbance  Statistical  Models 

The  statistics  of  the  disturbance  field  for  a  particu¬ 
lar  gravity  model  can  be  summarized  in  one  scalar  correla¬ 
tion  function  (see  Appendix  A).  With  this  function  as  a 


basis,  all  other  auto-  and  cross-correlations  of  anomalous 
gravity  terms  can  be  derived  CRef  43*  Field  theory  provides 
the  functional  interrelationships  between  these  terms 
CRef  3]-  Models  for  these  basic  functions  have  been  devel¬ 
oped  —  primarily  using  gravity  anomaly  autocorrelation 
function  as  a  basis.  Three  of  the  fully  developed  models 
are* 

a.  Linear  state  space  based  on  a  Gauss-Markov  process, 

b.  Anomaly  degree  variance  based  on  a  spherical 
harmonic  expression  of  the  correlation  function, 
and 

c.  Attenuated  white  noise  based  on  a  subterranean 
white  noise  process  for  anomalous  potential. 

Details  of  these  models  will  be  presented  in  examples 
of  Sections  IV,  V,  and  VI.  For  now,  the  point  is  that 
adequate  models  exist  to  support  the  statistical  evalua¬ 
tion. 


1.9  Present  Methods  of  Statistical  Analysis 

The  central  element  (Figure  2)  of  the  general  strategy 
is  the  analysis  method.  This  method  combines  the  three 
models  just  discussed  to  produce  an  estimate  of  the  navi¬ 
gation  system  error  statistics  caused  by  the  errors  in 
modeling  gravity.  Two  distinct  approaches  have  been  taken 
in  the  past:  Monte  Carlo  and  linear  state  space  covariance 
analysis.  These  methods  differ  considerably  and  must  be 
explored  separately  to  appreciate  the  inherent  strengths 
and  weaknesses  of  each. 
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1.9*1  Monte  Carlo  Method 

First,  the  Monte  Carlo  method  takes  a  direct  mission 
simulation  approach  (Figure  3).  It  is,  in  fact,  an  ensemble 
of  deterministic  cases  corresponding  to  an  ensemble  of 
residual  field  realizations.  The  gravity  disturbance 
profiles  are  produced  in  such  a  way  that  in  the  limit 
they  match  the  correlations  given  by  the  residual  field 
correlation  function.  These  simulations  may  be  either 
of  the  navigation  equations  directly  (whole-valued  simula¬ 
tion)  or  of  navigation  error  propagation.  In  either  case, 
the  ultimate  output  of  these  simulations  is  an  ensemble 
of  error- time  histories  from  which  means,  standard  devia¬ 
tions,  covariances,  and  other  sample  statistics  can  be  cal¬ 
culated.  Since  histograms  can  be  formed  for  any  particular 
error  at  any  point  in  time,  a  frequency  function  can  be 
approximated.  This  feature  is  not  of  significant  impor¬ 
tance  to  this  work  since  second  order  statistics  are  usually 
sufficient  to  specify  system  accuracy. 

The  significant  cost  of  the  Monte  Carlo  method  is  the 
computation  time  required  to  simulate  the  trajectories. 
Enough  cases  must  be  run  to  yield  statistics  of  high  confi¬ 
dence.  In  a  recent  advanced  cruise  missile  study,  Chatfield 
[Ref  10]  used  90  simulations  to  produce  the  desired  circu¬ 
lar-error-probable  statistic. 


While  this  cost  is  considerable,  the  fact  that  the 
analysis  could  be  performed  demonstrates  the  flexibility 
of  this  method.  The  Monte  Carlo  restrictions  on  trajec¬ 
tory  are  minor,  being  essentially  limited  to  the  ones 
imposed  by  the  trajectory  simulator.  In  addition  to  this 
scenario  flexibility,  considerable  freedom  exists  in  the 
functional  form  of  the  correlation  model.  The  gravity 
disturbance  profiles  are  produced  using  the  residual  field 
correlation  function.  While  the  residual  field  is  depen¬ 
dent  on  the  gravity  model,  there  are,  as  noted  above, 
several  developed  functional  models  to  summarize  the  corre¬ 
lations  of  the  residual  field.  The  Monte  Carlo  method 
imposes  no  limitation  on  the  functional  form  of  the  corre¬ 
lation  model. 

1.9.2  Linear  State  Space  Covariance  Analysis 

The  linear  state  space  covariance  approach  was  formu¬ 
lated  by  Levine  and  Gelb  in  a  landmark  paper  [Ref  11]  in 
1968.  This  method  and  the  attendant  stochastic  gravity 
model  have  been  refined  to  be  consistent  with  most  of  the 
field  theory  impositions  [Ref  12]. 

The  gravity  disturbance  terms,  u,  are  modeled  as  the 
output  of  a  linear  filter  driven  by  white,  gaussian,  inde¬ 
pendent  noise  sources.  The  state  vector  model  of  this 
linear  filter  is 

xg(t)  =  Fg('t)2ig('t)  +  Gg(t)q(t)  (4) 
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where  q(t)  is  a  white  gaussian  noise  process. 

The  output  equation  to  form  the  disturbance  terms  is 
u(t)  =  Cg(t)xg(t)  (5) 

The  white  noise  terms,  q( t) ,  are  modeled  to  have  correla¬ 
tion 

<fCq(t)qT(p)]  =  Qg(  t)  6(p-t)  (6) 


where  6(-)  is  the  Dirac  delta  function.  An  augmented 
state  (see  Figure  4)  can  now  be  defined  as 


The  augmented  state  vector  differential  equation  is 

X&(  t)  =  Fa(t)xa(t)  +  Ga(t)q(t)  (8) 

where 


patt) 


and 


F(t)  . 

G(  t)  Cg(  t)  * 

0 

Fg(t)  _ 

(9) 


Ga(t) 


0 


°g<t) 


(10) 


With  this  augmented  form  the  system  covariance  integral 


becomes 

Pa(t)  =  St  St  $a(t,p)Ga(p)Qg(p)  6(q-p)Ga(q)^a(  t,q)dpdq 

0  0  (11) 

where 

=  Fa(t)$a(t,t1)  (12a) 


1? 


LINEAR  STATE  SPACE  COVARIANCE  ANALYSIS 


Figure  4.  Linear  State  Space  Covariance  Method 
Model  Interface 
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subject  to 

®a(W  =  1  (12b) 

The  delta  function  allows  a  reduction  to 

?a( ‘t)  =  $a(t,p)Ga(p)Q  (p)G^(p)^(t,p)dp  (13) 

Applying  Leibnitz’  rule  now  yields 

Pa(t)  =  Fa(t)Pa(t)+Pa(t)FT(t)  +  Ga(t)Qg(t)Ga(t)  (14) 

This  equation  can  be  solved  numerically  by  a  number  of 
powerful  and  efficient  numerical  integration  techniques. 

The  fact  that  the  solution  can  be  cast  in  this  mold  is 
the  most  persuasive  argument  for  using  linear  state  space 
methods . 

The  key  to  this  efficiency  is  the  delta  function 
which  allows  one  level  of  integration  to  be  solved  directly 
without  approximation.  The  Gauss-Markov  model  which  pro¬ 
duced  this  advantage  is  the  root  of  errors  that  can  result 
when  linear  state  space  covariance  analysis  is  used  on 
complex  trajectories. 

Obviously,  the  linear  state  space  covariance  analysis 
restricts  the  gravity  disturbance  statistical  model  form  to 
one  compatible  with  linear  state  space  description.  This 
restriction,  while  undesirable,  is  much  less  serious  than 
the  loss  of  trajectory  generality  which  accompanies  the 
time-domain  use  of  the  Gauss-Markov  model.  Gravity  disturb¬ 
ance  is  a  spatial  function  and  hence  the  statistical  repre¬ 
sentation  is  a  spatial  process.  The  conversion  from 
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spatial  to  a  temporal  representation  of  the  statistical 
process,  symbolized  in  Figure  5i  causes  the  resulting 
model  to  err  when  the  trajectory  varies  from  a  particular 
form. 

The  details  of  the  spatial  to  temporal  model  conver¬ 
sion  are  presented  in  Appendix  B.  To  understand  why  this 
conversion  creates  error,  one  need  only  consider  the  vari¬ 
ables  involved.  The  spatial  model  has  central  angle,  \Jj, 
as  an  independent  variable.  Central  angle  is  not  a  scalar 
since,  in  general,  the  total  central  angle  change  is  not 
equal  to  the  sum  of  the  central  angle  changes  along  mission 
segments.  Time  is  a  scalar  and  the  conversion  to  the  time 
independent  variable  means  that  the  underlying  relationship 
is  mismodeled,  in  general.  For  the  special  case  of  a  great 
circle  trajectory,  the  total  central  angle  change  is  the 
sum  of  the  segment  changes  and  the  conversion  is  also  only 
strictly  valid  for  great  circle  missions  with  non-decreasing 
central  angle  of  less  than  180°. 

Another  trajectory  restriction  is  imposed  by  the 
original  spatial  Gauss-Markov  model.  This  formulation  does 
not  treat  changes  in  altitude,  and  therefore,  the  model 
is  only  valid  for  constant  altitude  trajectories  as  well. 
Equivalent  statistics  may  be  calculated  from  the  original 
at  other  altitude  levels  [Ref  133;  however,  combining  these 
into  a  linear  system  format  requires  some  improvisation. 
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In  summary  then,  the  only  trajectory  for  which  the 
linear  state  space  covariance  analysis  is  valid  is  a  con¬ 
stant  altitude,  great  circle  case  with  a  monotonically 
non-decreasing  central  angle  of  less  than  180°.  So,  the 
considerable  computation  efficiency  achieved  through  (14) 
has  the  severe  trajectory  limitations  as  a  price.  Although 
the  trajectory  restrictions  of  the  Gauss-Markov  model 
are  well-known,  the  extent  of  errors  induced  when  the 
trajectory  restrictions  are  violated  is  not  known. 

Numerical  examples  will  be  presented  in  Section  IV  to 
explore  this  area. 

1.10  Motivation  for  Research 

Adequate  models  exist  for  trajectory,  for  navigation 
system,  and  for  gravity  disturbance  statistics,  but  the  two 
present  methods  of  statistical  analysis  have  serious  dis¬ 
advantages  which  limit  their  usefulness.  The  Monte  Carlo 
method  has  virtually  no  restrictions,  but  the  requirement 
to  produce  samples  until  the  statistics  stabilize  is  often 
costly.  On  the  other  hand,  the  trajectory  restrictions 
of  the  linear  state  space  covariance  method  significantly 
limit  the  usefulness  of  this  method.  The  gravity  disturb¬ 
ance  statistical  model  is  also  limited  to  a  particular 
Gauss-Markov  form  for  the  linear  state  space  method,  and 
this  constraint  might  limit  the  utility  of  this  evaluation 
method  on  some  alternative  on-line  gravity  models  such  as 
a  truncated  spherical  harmonic  model. 
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An  alternative  to  these  methods  is  needed  which  is 
more  flexible  than  linear  state  space  covariance  analysis 
and  which  is  more  efficient  than  Monte  Carlo.  To  develop 
such  an  alternative  is  the  objective  of  this  research. 

1.11  Overview  of  the  Research 

The  following  sections  document  the  research  to 
develop,  verify  and  demonstrate  the  alternative  statistical 
analysis  method.  Section  II  presents  the  theoretical 
development  and  Section  III  the  numerical  method.  The  new 
method  is  compared  to  the  linear  state  space  method  in 
Section  IV  and  to  Monte  Carlo  in  Section  V.  These  compari¬ 
sons  constitute  a  thorough  verification  of  both  the  theo¬ 
retical  and  numerical  aspects  of  the  new  analysis  technique. 
Demonstrations  are  then  presented  (Sections  VI,  VII,  and 

VIII)  which  show  the  flexibility  of  this  new  analysis  method 
to  treat  a  variety  of  gravity  models  and  gravity  disturb¬ 
ance  statistical  modeling  methods.  In  a  final  develop¬ 
ment,  the  effects  of  Kalman  filter  updates  on  gravity 
induced  navigation  error  statistics  are  included  (Section 

IX)  in  the  analysis  algorithm,  and  this  update  process  is 
verified.  In  all,  a  new  statistical  analysis  technique 
is  presented  to  evaluate  the  system  errors  induced  by 
gravity  modeling  errors. 
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II.  Development  of  a  Gravity  Model  Performance 

Measure  Based  on  an  Inertial  Navigation  Error 
Covariance  Integral _ 

In  Section  I,  the  present  methods  of  statistical 
analysis  were  reviewed.  The  Monte  Carlo  method,  while 
flexible,  has  significant  computational  cost.  The  linear 
state  space  covariance  analysis  is  computationally  efficient 
but  has  unacceptable  trajectory  restrictions  associated 
with  the  Gauss-Markov  gravity  disturbance  statistical 
model.  To  develop  an  efficient  yet  flexible  alternative 
analysis  method,  a  line  of  development  is  suggested  which 
is  parallel  to  the  linear  state  space  (Figure  5)  but  which 
avoids  the  dependence  on  the  Gauss-Markov  statistical  model. 
The  purpose  of  this  section  is  to  present  the  theoretical 
development  of  the  covariance  integral  as  a  candidate 
alternative  statistical  analysis  method. 

2.1  General  Approach 

The  genealogy  of  the  linear  state  space  covariance 
analysis  method  (portrayed  symbolically  in  Figure  5)  and 
the  parallel  approach  to  be  taken  here  are  shown  in  Figure 
6.  The  crucial  change-of- variables  (from  central  angle  ifr 
to  time  t)  does  not  occur  in  this  new  "covariance  integral" 
formulations  the  spatial  dependence  of  the  correlations  is 
retained.  This  approach  does  not  exclude  the  Gauss-Markov 
gravity  disturbance  statistical  model  from  use.  But,  it 
does  mean  that  the  gravity  disturbance  correlations  from 
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that  model  are  used  directly,  and  not  the  temporal  linear 
filter  model  of  gravity  disturbance  (4).  As  a  consequence, 
the  integral  expression  of  the  navigation  error  covariance 
will  not  contain  a  simplifying  Dirac  delta  function  as  in 
(11).  The  covariance  integral  becomes,  then,  the  end 
object  of  the  theoretical  development.  The  final  analysis 
will  be  simply  a  numerical  approximation  of  this  integral. 
Theoretical  development  of  the  new  "covariance  integrals" 
statistical  evaluation  method  is  made  in  the  next  subsec¬ 
tion. 

2. 2  The  Covariance  Integral  Theoretical  Development 

The  sample  space  definition  is  the  key  to  any  statis¬ 
tical  analysis.  The  view  taken  here  is  that  the  desired 
error  statistics  represent,  in  a  performance  index  sense, 
the  expected  performance  over  some  range  0  of  possible 
missions.  Then  to  produce  statistics  representative  of 
missions  in  0,  the  statistical  expectation  must  be  over 
all  06 0 .  A  mission  0  is  simply  a  position-time  history  or 
trajectory  and  corresponds  to  one  sample  from  the  space  0. 
For  each  0,  a  navigation  propagation  characteristic  is 
defined  and  a  gravity  disturbance  is  generated.  This 
disturbance  is  implicitly  a  function  of  time  being  expli¬ 
citly  dependent  on  the  r(t)  dictated  by  0.  Then  the  navi¬ 
gation  error  dynamics  (3)  may  be  written  as 

x(t;0)  =  F( t ; 9)  x(t|0)  +  G(t;0)  u(t»0)  (15) 
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with  initial  error  state  x(tQ}0)  which  will  be  retained 
in  the  development  for  generality. 


Since  (15)  is  a  linear  system,  there  exists  a  state 
transition  matrix  $(t,t^;9)  which  satisfies 

SKt.t^e)  =  f( t ; 0)  ^(t.t^e)  (16) 

and 

$(t,tj0)  =  I  (17) 

This  matrix  also  satisfies  the  semigroup  property 

^(t^.t-j^e)  =  $(t^,t25  0)  ^(tg.^ie)  (18) 

which  will  be  essential  in  the  upcoming  numerical  analysis. 

For  now  these  properties  allow  the  solution  of  (15)  to  be 

cast  in  the  form 
t 

x(t}0)  =  ft  $>(t,p;0)  G(  p ; 0)  u(p;0)  dp  +  $>(  t,  tQ ;  0)  x(  tQ ;  0) 

°  (19) 

The  next  step  toward  forming  the  covariance  matrix  is  to 
produce  the  outer  product 

x(t50)xT(t;0)  =  <2>(t,p;0)G(p;0)  [u(p; 0)uT(q; 0) ] 

o  o 

GT( q  j 0)$T( t,q; 0) dpdq 


+  {$(t,to50)  [x(to;0)uT(p;0)]  GT(  p ;  0) 

$T(t,pj  0)dp| 

+  {©(t,to»0)  s\  [x(tos0)uT(p;0)]  GT(p;0) 
$  T(t,p;0)dp|T 

+  $(t,to»0)  Cx(to;0)xT(to»0)]  $T(t,tos0) 


The  navigation  error  state  covariance  P(t)  is  defined  as 
the  expectation  of  this  outer  product  over  all  mission  0 
within  the  mission  region  0.  That  is, 

P(t)  =  $  fx(t;0)xT(t;0)]  (21) 

0  e  e  L  J 

Now  to  cast  the  problem  in  a  form  for  which  computer 
solution  is  practical,  assume  that  the  mission  region  0 
has  been  selected  sufficiently  small  that  one  mission,  the 
design  mission  0Q,  can  represent  the  mission  geometries  and 
the  navigation  error  dynamics  throughout  that  region  (Figure 
7).  The  disturbance,  u(t,0),  varies  more  with  geographic 
variations;  hence,  is  retained  in  that  form.  The  design 
mission  9o  defines  an  occurrence  of  a  position- time  history 
r(t)  for  t€(tQ,tf).  The  present  assumptions  mean  the 
following  approximation  replacements  can  occur  as  defined 
below 

F(t;0)^-F(t)  =  F(t;0o)  (22) 

G(tse)«*-G(t)  =  G(t;60)  (23) 

$(  t,  tx;  0)-*-*$(t,  t^  =  $(t,t1;0o)  (24) 

This  replacement  will  simplify  taking  the  expectation 
since  F( t) ,  G(t),  and  $(t,p)  are  not  changed  over  the 
range  of  0.  Without  such  a  simplification,  the  expectation 
expressed  by  (21)  would  be  extremely  difficult  to  approxi¬ 
mate  by  other  than  Monte  Carlo  analysis.  Applying  (21)  to 
(20)  then  yields 
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Figure  7-  Covariance  Integral  Mission  Sample 
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P(t)  =  /*  /+  $>(t,p)G(p)  $  [u(pje)uT(q»e)]  GT(q)$T(t,q)dpdq 

**o  o  a«e 


+  E(t)  +  ET(t)  +  $(t,to)P(to)0T(t,to)  (25) 


and 

E(t) 


=  $(t,t0)J^  $  [x('to!0)iiT(p;e)'|  GT(p)$T( t,p)dp 

o  ete  *-  J 


(26) 


The  initial  covariance,  P(tQ),  and  hence,  x(tQ;e) ,  will 
he  modeled  as  identically  zero  for  the  following  develop¬ 
ment.  A  perfect  initial  condition  is  thus  assumed.  The 
penalty  for  this  simplification  is  that  calibration  and 
alignment  errors  which  are  geodetic  in  origin  cannot  be 
included  as  initial  conditions.  Should  it  become  important 
to  include  such  terms  with  the  gravity  modeling  error  evalu¬ 
ation,  some  additional  development  will  be  required  to 
approximate  (26) .  Alternatively,  calibration  and  align¬ 
ment  can  be  included  as  part  of  the  design  mission.  The 
navigation  error  propagation  model  and  the  trajectory  must 
reflect  this  fact. 

By  the  above  simplifying  assumptions,  the  last  three 
terms  in  (25)  are  removed  yielding 

P(t)  =  ft  S \  $(  t,p)G(p)Q(p,q)GT(q)$T(  t,q)dpdq  (2?) 

o  '’o 

where 

Q(p,q)  =  (u(pj  9)  uT(q;  0)  }  (28) 

6(0 

This  geodetic  error  correlation  function  Q(p,q)  sum¬ 
marizes  the  statistical  relations  between  gravity  errors 
at  time  p  with  those  at  time  q.  Since  u(p,6)  is,  in 


reality,  u[r(p)]  for  r(p)  defined  by  9,  it  is  more  correct 
to  express  this  function  in  terms  of  positions  (i.e., 

QCr(p) , r(q) ]) •  Again,  0Q  yields  the  required  position 
history,  and  r(p)  and  r(q)  are  defined  by  0Q  for  evalua¬ 
tions  of  the  correlation  function  model. 

The  elements  of  the  Q-matrix  are  auto-  and  cross¬ 
correlation  functions  which,  as  previously  discussed,  can 
be  derived  from  one  scalar  correlation  function  model  £see 
Appendix  A].  For  this  work  it  suffices  to  note  that  a 
body  of  theory  exists  on  the  relationships  of  the  corre¬ 
lation  of  anomalous  gravity  disturbance  vector  6g  and 
undulation  of  the  geoid  N.  Related  terms  are  vertical 
deflections,  ij  and  f  ,  gravity  anomaly  Ag,  and  anomalous 
potential  T.  These  anomalous  gravity  terms  are  defined  in 
Appendix  A.  One  purpose  here  is  to  develop  a  scheme  gene¬ 
ral  enough  to  allow  a  correlation  modeling  choice. 

Equation  (2?)  is  the  desired  expression  of  navigation 
error  covariance  as  an  integral  involving  the  trajectory, 
error  propagation,  and  correlation  function  models.  Numeri¬ 
cal  methods  will  be  derived  in  the  next  section  to  approxi¬ 
mate  (27)  directly.  Another  formulation  is  suggested  by 
the  linear  state  space  methods. 

Equation  (27)  is  similar  to  (11)  of  the  linear  state 
space  statistical  analyses.  Computational  advantage  is 
gained  in  that  development  when  Leibnitz*  rule  is  applied. 
Following  this  example  yields 
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P(t)  =  F(t)P(t)  +  P(t)FT(t)  +  G(t)D(t)  +  DT(t)GT(t)  (29a) 


where 

D(t)  =  /!  Q(t,p)GT(p)  $T(t,p)dp  (29b) 

xo 

Since  time  appears  as  an  argument  of  Q( • , • )  differentia¬ 
tion  of  (29b)  will  not  yield  a  simple  form  in  general. 
Combining  (29b)  with  the  integral  form  of  (29a)  yields  a 
pair  of  nested  single  integrals  which  can  be  approximated, 
instead  of  directly  approximating  (27) •  Approximations 
based  on  this  nested  integrals  approach  will  also  be 
developed  in  the  next  section. 

2.3  Comparisons  to  Present  Statistical  Analysis  Methods 
The  objectives  of  this  "covariance  integral"  evalua¬ 
tion  method  are  efficiency  and  flexibility  advantage  over 
existing  methods  of  gravity  model  statistical  evaluation. 
The  efficiency  of  a  numerically  approximated  covariance 
integral  is  compared  to  Monte  Carlo  in  Section  V.  This 
comparison  must  wait  until  then  because  first  the  develop¬ 
ment  of  numerical  methods  for  approximating  the  covariance 
integral  must  be  accomplished.  These  numerical  methods  are 
developed  in  Section  III.  However,  the  flexibility  compari 
sons  can  be  made  at  this  point. 

Flexibility,  as  used  here,  is  the  facility  to  process 
a  variety  of  models  in  either  of  the  three  categories: 
trajectory,  navigation  system,  and  disturbance  statistics. 


The  model  flexibility  for  the  new  covariance  integral 
formulation  is  presented  in  Table  I  along  with  the  two 
current  analysis  methods.  This  new  method  is,  as  desired, 
more  flexible  than  the  linear  state  space  covariance 
analysis. 

The  trajectory  flexibility  is  the  major  concern  with 
the  linear  state  space  method.  The  formulation  of  (27), 
or  the  equivalent  (29),  impose  no  special  constraint  on 
r(t),  the  trajectory.  As  long  as  the  Q-matrix  and  F-matrix 
can  be  evaluated  using  r(t),  the  covariance  integral  offers 
a  means  of  evaluating  the  system  performance  of  the  gravity 
model  on  any  trajectory. 

The  linear  state  space  evaluation  approach  is  forced 
to  use  the  Gauss-Markov  gravity  disturbance  statistical 
model.  This  restriction  is  not,  of  itself,  a  major  concern, 
but  it  is  a  definite  limitation.  On  the  other  hand,  the 
new  formulation  requires  only  that  the  disturbance  correla¬ 
tions  be  modeled  as  spatially  dependent  quantities  for 
which  r(p)  and  r(q)  from  the  design  mission  will  serve  as 
inputs.  This  flexibility  allows  the  analyst  options  in  this 
modeling  which  may  simplify  the  overall  task.  A  demonstra¬ 
tion  using  three  different  model  forms  for  the  gravity  dis¬ 
turbance  statistics  is  presented  in  Section  VI. 

The  navigation  model  is  restricted  in  this  new  formula¬ 
tion.  The  linear  state  space  navigation  error  propagation 
model  shared  by  linear  state  space  and  the  covariance  integral 
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formulations  is  not  considered  a  serious  drawback  because 
the  model  is  widely-accepted  as  a  valid  system  description. 

The  covariance  integral  approach  gives  the  desired 
model  flexibility.  The  promise  of  greater  efficiency 
than  Monte  Carlo  can  be  shown  only  after  a  numerical 
approximation  algorithm  is  developed  and  proven. 


III.  Numerical  Approximation  the  Covariance  Integral 

In  Section  I,  the  motivation  was  given  for  developing 
a  new  statistical  method  to  evaluate  system  accuracy  per¬ 
formance  of  gravity  models.  The  theoretical  basis  for  the 
covariance  integral  (27).  presented  in  Section  II,  offers 
considerable  freedom  in  model  selections.  Except  for  path¬ 
ologically  simple  cases,  a  closed-form  expression  for  covari¬ 
ance  integral  cannot  be  expected,  so  a  means  of  approxi¬ 
mating  the  integrations  is  needed.  The  main  purpose  of 
this  section  is  to  investigate  alternative  approximation 
procedures  and  to  select  one  procedure  for  futher  verifica¬ 
tion,  demonstration,  and  development.  Approximating  the 
covariance  integral  is  straightforward,  but  computational 
cost  becomes  prohibitive  for  realistic  problems.  Therefore, 
this  section  also  addresses  the  data  logistics  issue  which 
is  critical  in  forming  an  economically  feasible  analysis 
method. 

Two  general  lines  can  be  followed  to  approximate  the 
covariance  integral.  A  direct  approximation  can  be  made 
by  replacing  the  double  integration  of  (27)  with  finite 
double  summations.  An  alternative,  referred  to  as  the 
Nested  Integrals  method,  is  to  approximate  the  integrals 
of  (29b)  and  (29a)  in  that  order.  Recursive  algorithms 
will  be  developed  along  both  these  theoretically  equiva¬ 
lent  lines  after  the  common  computational  problems  are 
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addressed.  Both  the  accuracy  and  efficiency  of  these 
candidate  algorithms  will  then  be  assessed  in  order  to 
select  one  method  for  further  study. 

The  next  three  subsections  are  devoted  to  the  computa¬ 
tional  details  of  these  algorithms.  The  reader  not 
interested  in  these  matters  may  wish  to  proceed  to  sub¬ 
section  3.4  where  the  comparison  of  the  developed  algo¬ 
rithms  begins. 

3.1  Common  Computational  Problems 

The  theoretical  work  of  this  research  could  be  covered 
without  discussing  the  techniques  used  in  producing  numeri¬ 
cal  results.  Such  treatment  would  neglect  an  issue  criti¬ 
cal  to  the  use  of  the  covariance  integral  statistical 
method.  In  developing  computer  programs  to  execute  the 
approximations  of  (27)  and  (29).  considerable  study  and 
experiment  went  into  the  data  flow  and  algorithmic  struc¬ 
ture.  This  subsection  presents  the  details  of  these  soft¬ 
ware  design  decisions,  because  the  solution  reached  in 
this  study  may  be  applicable  to  other  studies. 

Several  arbitrary  decisions  must  be  made  in  designing 
software  for  a  general  requirement  like  "Evaluate  (27)  with 
a  numerical  approximation  over  a  wide  range  of  trajectory, 
error  propagation,  and  gravity  disturbance  models." 

Other  general  specifications  include  efficiency  and  accu¬ 
racy  requirements.  Factors  which  influence  these  decisions 
are  1 
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a.  The  specific  problem  to  be  solved.  The  strategy 
used  on  simple  problems  like  the  Schuler  loop  examples, 
which  are  presented  later  in  this  section,  is  certainly 
different  from  that  used  on  the  full-scale  problems  pre¬ 
sented  in  the  following  sections. 

b.  Computing  system  capacities.  Available  computer 
core,  input-output  efficiency,  and  throughput  are  examples 
of  the  capacities  of  a  specific  computational  facility 
which  should  be  considered. 

c.  Available  support  software.  Computer  file  record 
management  and  subroutine  libraries  should  be  exploited. 

d.  General  user  environment.  Priority  is  established 
at  most  large  computing  facilities  based  on  many  factors, 
and  this  general  environment  influences  the  nature  of  the 
software  design. 

The  decisions  made  in  preparing  software  for  this  study 
were  based  on  using  a  CDC  CYBER  7^  computing  facility  in 
an  environment  which  encourages  short  multi-step  processes 
over  long  computation  time  jobs.  The  following  decisions 
were  influenced  by  this  environment  and  this  inherent 
bias  should  be  remembered.  With  this  warning,  the  following 
discussion  gives  one  example  of  a  solution  to  the  computa¬ 
tional  burden  associated  with  approximating  the  covariance 
integral. 

Either  a  direct  approximation  of  (27)  or  the  Nested 
Integrals  approximation  of  (29)  requires  r(t)  and  <5(p,q) 
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data.  Methods  for  producing  r(t)  exist  CRef  7  ]  for  trajec 
tories  too  complex  for  closed-form  mathematical  expressions 
The  ^-matrices  can  be  produced  using  (16)  and  (17)  by  ordi¬ 
nary  differential  equation,  predictor-corrector  methods 
Ce.g.  Ref  14].  Once  produced,  the  logistics  of  supplying 
these  data  must  be  resolved.  This  issue  is  treated  first 
because  the  algorithmic  forms  are  influenced  by  the  data 
flow.  Since  the  Nested  Integrals  data  requirements  are 
practically  the  same  as  the  direct  approximations,  the 
direct  method  will  be  used  as  an  example  in  the  following 
discussion. 

The  direct  integral  approximations  will  be  based  on 
a  finite  double  summation  replacing  (27).  That  is, 

P(n)  =  2  Z  $  (n,i)G(i)Q(i,  j)GT(  ,j)$T(n,  j)  At.At.s  (i,j) 

i=0  j=0  1  J  n 

(30) 

where  s  (i,j)  is  the  quadrature  weighting  associated  with 
the  integrand  evaluation  at  (p,q)  =  (i  ,  j  ).  Details  of 
the  weighting  factor  are  covered  later  in  the  algorithm 
developments.  Attention  will  now  be  given  to  the  nature 
of  the  data  and  to  the  manner  in  which  it  is  produced. 

Fundamental  to  the  formation  of  (30)  is  a  discreti¬ 
zation  of  the  time  interval  (tQftn).  One  can  view  the 
integrals  of  (27)  as  a  process  applied  to  a  signal,  the 
integrand.  With  this  approach,  the  time  steps  should  be 
small  enough  to  give  at  least  adequate  representation  of 
the  integrand  as  specified  by  Shannon’s  sampling  theorem. 


The  frequency  spectrum  of  the  integrand  is  affected  by 
all  three  analysis  models.  The  trajectory  affects  both 
the  $  and  Q  matrices.  The  inertial  navigation  error  propa¬ 
gation  model  along  with  the  design  mission  trajectory  model 
define  the  F  matrix  used  in  computing  $.  The  gravity 
disturbance  correlation  model  is  embodied  in  Q,  which  has 
position  arguments  defined  by  the  trajectory.  The  dynamics 
of  the  trajectory,  the  dynamics  of  the  error  propagation 
model,  and  the  spatial  spectrum  of  the  gravity  disturbance 
correlation  model  must  all  be  considered  in  subdividing 
the  time  interval. 

Once  a  proper  time  partition  is  defined,  the  integrand 
must  be  evaluated  for  every  pair  (t.  ,t.)  such  that  both  t. 

•1.  J  -L 

and  tj  are  in  the  set  {t0,t1 . tN).  This  $ -matrix 

data  need  is  a  critical  computational  issue  which  will  be 
addressed  at  length. 

If  n  is  the  number  of  error  states,  the  integral  sum- 
s 

2 

mation  can  also  be  reduced  from  n  terms  to  £n  (n  +1) 

s  s  s 

because  P(t)  is  symmetric.  This  symmetric  property,  as 
will  be  shown  later,  can  be  used  to  decrease  the  number  of 
$(i,j)  matrices.  That  is,  the  only  $(i,j)  matrices  required 
are  for  (i  ,  j  )  pairs  such  that  0  <  j  <  i  <  N.  Since 
$(i,j)  =  I  when  i=j,  the  true  requirement  is  for  (  i,  j) 
pairs  such  that  0  <  j  <  i  <  N.  For  example,  $(4,7) 
need  not  be  calculated,  whereas  $(7,4)  is  required.  The 
number  of  $-matrices  required  is  £(N+l)(N+2)  -  (N+l)  - 
|N(N+1) .  To  demonstrate  the  magnitude  of  this  data  burden, 
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consider  a  500-point  time  partition  and  9X9  $-ma trices. 

The  12,250  $-ma trices  would  require  10,145,250  words  of 
storage. 

Clearly,  the  logistics  of  these  data  must  he  care¬ 
fully  planned  to  make  this  analysis  method  economically 
feasible.  The  semi-group  property,  equation  ( 18),  can  be 
used  to  reduce  this  data  storage,  but  at  a  cost  in  compu¬ 
tation  time.  For  discussion  purposes,  let  N=3-  The  approxi¬ 
mation  stated  in  (30)  requires  $(1,0),  $(2,1),  $(2,0), 

$(3>2),  $(3»1),  and  $(3,0).  All  of  these  can  be  produced 
from  $(3,2),  $(2,1),  and  $(1,0)  using  (18).  To  see  this, 
note  that 


$(1,0)  =  $(1,0) 

(31a) 

$(2.1)  =  $(2,1) 

(31b) 

$(2,0)  =  $(2,1) 

$(1,0) 

(31c) 

$(3,2)  =  $(3,2) 

(3ld) 

$(3,1)  =  $(3,2) 

$(2,1) 

(3le) 

$(3,0)  =  $(3,2) 

$(2,1)  $(1,0) 

(3lf) 

For  N=3,  then,  the  storage  requirement  can  be  reduced  from 
six  to  three  $-ma trices.  In  general,  this  technique 
requires  only  N  $-matrices  as  opposed  to  the  sN(N+l) 
needed  for  (30) .  For  the  500-point  example  with  9x9 
$  -matrices,  the  storage  requirement  decreases  from  10,145, 
250  to  40,500  words. 
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Note  that  the  production  of  $(3,0)  above,  requires 
two  matrix  multiplications.  This  computational  cost  is 
offset  by  the  fact  that  $(3,0)  is  not  produced  by  solving 
(16)  initialized  with  $(0,0)  =  I.  Note  that  although  the 
product  sequence  runs  backward  in  time,  each  $(i+l,i) 
results  from  a  forward-time  integration  from  the  initial 
condition  $(i,i)  =  I.  Questions  of  numerical  stability 
must  be  asked  concerning  the  possibly  hundreds  of  matrix 
multiplications  which  full-scale  problems  might  require. 

For  the  studies  presented  in  this  and  later  sections,  a 
simple  test  was  performed  to  test  this  $-matrix  production 
technique.  $(N,0)  was  produced  by  two  methods.  First, 
the  $(i,i-l)  matrices  were  calculated  using  (16)  and  (17). 
From  these, 

$(N,0)  =  $(N,N-1)  $(N,N-2)  ...  $(2,1)  $(1,0)  (32) 

was  calculated.  Then,  $(N,0)  was  produced  by  integrating 
(16)  initialized  with  $(0,0)  =  I.  The  elements  of  these 
two  alternative  $(N,0)  matrices  were  individually  compared. 
The  differences  observed  were  on  the  order  of  the  computer 
word  accuracy  limit  and  this  technique  was  judged  accurate 
enough  for  the  studies  which  follow. 

These  test  results  are  unique  to  the  computing  system 
and  model  choices.  For  other  applications,  a  test  similar 
to  the  one  above  is  advised  prior  to  conducting  an  analysis 
using  the  technique  of  repeated  multiplications. 


It  should  be  clear  from  ( 3 1)  that  a  preferred  order 
of  computation  exists.  Repetitious  calculations  can  be 
avoided,  for  example,  if  $(3»1)  is  calculated  and  used 
prior  to  $(3,0).  The  associative  and  commutative  proper¬ 
ties  of  matrix  addition  permit  the  summation  of  (30) 
such  that  $(3,2)  comes  before  $(3,1)  which  corner;  before 
$  (3,0)*  In  general,  an  algorithmic  form  is  desired  in 
which  $(n,n-l) ,  $(n,n-2) ,  $(n,0)  are  required  in  that 

order.  The  second  index  decreases  —  moves  backward  in 
time.  This  order  of  summation  is  therefore  backward  from 
the  way  the  sum  would  normally  be  formed.  Explicit  refer¬ 
ence  to  this  reverse-chronological  order  will  be  made  at 
those  points  in  the  algorithm  developments  where  the  sum¬ 
mations  are  modified  to  this  desired  order. 

Another  data  logistics  issue  is  whether  to  produce  the 
$ -matrices  in  parallel  with  P  computations  or  to  produce 
all  the  data  before  computing  any  P  values.  The  parallel 
production  method  means  $(n,n-l)  would  be  produced  and 
stored  just  prior  to  the  P(n)  computation.  All  previous 
$ -matrices  would  have  been  produced  on  previous  steps  and 
stored  for  retrieval.  This  approach,  therefore,  almost 
forces  a  chronologically  sequential  $-matrix  storage.  This 
storage  mode  would  be  prohibitively  costly  in  computer 
input-output  time  if  the  reverse-chronological  summation 
form  is  used.  Another  disadvantage  is  that  the  trajectory 
model  which  supports  $  computation  would  either  have  to 


reside  in  computer  core  with  the  P  calculation  or  be 
overlayed  iteratively  with  the  P  calculation.  Except  for 
simple  analyses,  the  computer  core  requirements  would  be  a 
severe  disadvantage  on  one  hand  or  the  peripheral  proces¬ 
sing  would  be  costly  on  the  other. 

If  the  parallel  computing  structure  has  disadvantages, 
the  serial  design  has  distinct  positive  attributes.  First, 
the  trajectory  generation  and  state  transition  matrix  cal¬ 
culation  form  a  natural  partition  for  the  analysis.  These 
two  data  sets  provide  all  the  data  for  (30)  except  that 
associated  with  the  correlation  model.  The  <l>-matrix  calcu¬ 
lation  requires  the  trajectory  model  already,  so  a  parallel 
calculation  imposes  no  additional  computer  core  costs.  The 
data  from  both  trajectory  and  $  calculations  can  be  gene¬ 
rated  in  a  natural  ascending  time  order,  and  the  data  order 
can  be  reversed  in  an  inexpensive  data  processing  step. 

Once  these  data  have  been  produced  and  filed,  they  can  be 
used  repetitively,  if  needed,  without  the  cost  of  regene¬ 
ration. 

Now,  attention  must  turn  to  the  data  retrieval  costs. 
For  reasonably  sized  problems  like  the  examples  in  later 
sections,  a  straightforward  approach  to  retrieving  these 
data  for  the  required  computation  could  require  several 
hours  of  peripheral  processing  time  due  to  the  input-output 
operations  involved.  This  cost  is  likely  to  be  the  dominant 
factor  in  designing  any  software  to  approximate  the  covari¬ 
ance  integral.  To  understand  the  magnitude  of  this 


prospective  cost,  consider  the  following.  To  compute 
P(l),  $(1,0)  is  needed;  to  produce  P(2),  $(2,1)  and  $(1,0) 
are  required;  and  to  form  P(N)  a  total  of  N  matrices  are 
retrieved.  In  computing  P(l),  P(2),  ...,  P(N-l) ,  and 
P(N) ,  a  total  of  |-N(N+1)  $-matrix  fetch  operations  must  be 
performed.  In  later  sections,  example  problems  are  treated 
for  which  N=4o8.  For  such  problems,  a  total  of  83,436 
fetch  operations  are  required.  Since  these  matrices  can¬ 
not  all  reside  in  computer  core,  the  input-output  time  can 
be  substantial. 

Methods  of  reducing  this  computational  burden  must  be 
considered.  Even  though  all  $-ma trices  cannot  be  in  core 
simultaneously,  certainly  several  may  be.  To  the  existing 
data  processing  step,  a  task  can  be  added  to  pack  the 
$-matrices.  Several  chronologically  contiguous  $-ma trices 
could  be  packed  into  each  record.  Then,  several  matrices 
would  be  brought  into  core  with  each  fetch  operation.  If 
k  matrices  are  packed  into  each  record  and  if  N/k  is  an 
integer,  the  number  of  record  retrievals  is  decreased  from 
£N(N+1)  to  |N(N/k+l) .  With  N=408  and  with  k=4,  the  83,436 
record  fetches  are  reduced  to  21,012.  Even  though  more 
data  is  transferred,  a  savings  of  input-output  time  is 
realized  from  the  avoided  fetch  operations. 

The  cost  in  core  storage  for  this  input-output  effi¬ 
ciency  is  modest.  Instead  of  storing  one  $-matrix  in 
core  at  a  time,  k  of  them  must  be  stored.  The  computer 
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core  storage  increases  from  n  to  kn  .  The  necessary 
trade-off  between  core  storage  and  input-output  costs  is 
a  subjective  decision  influenced  primarily  by  the  computing 
environment  priority  system.  The  method  of  packing  records 
offers  one  means  of  lowering  input-output  time  should  such 
economy  be  desirable.  If  the  experience  of  these  studies 
is  a  general  indicator,  this  issue  will  be  crucial  in 
using  a  covariance  integral  statistical  analysis. 

The  records  used  in  these  studies  were  placed  in 
reverse  chronological  order,  and  four  nascent  records  were 
packed  into  one  record  prior  to  the  filing  step.  The  type 
of  file  structure  has  not  been  addressed  in  the  previous 
discussions.  The  filing  suggested  is  sequential,  but  a 
simple  sequential  file  structure  was  not  used.  Three  file 
structures  were  considered:  sequential,  random  access  or 
indexed,  and  indexed-sequential.  The  decision  to  use  the 
indexed-sequential  structure  was  based  partly  on  experience 
and  partly  on  experiment.  To  clarify  the  rationale  for  this 
choice,  these  file  structures  will  be  discussed  in  terms 
of  how  the  integral  approximation  would  proceed  using  each 
different  structure. 

a.  Sequential  file.  This  file  structure  is  probably 
the  most  familiar  since  the  normal  FORTRAN  read  and  write 
operations  use  and  create  sequential  files.  Record  retrie¬ 
vals  are  sequential  and  usually  proceed  one  record  at  a 
time.  For  the  application  here,  the  (^-matrices  would  be 
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stored  in  reverse  chronological  order  —  each  record  being 
packed  with  k  data  sets.  To  calculate  P(n) ,  the  first 
record  retrieved  should  contain  the  $(n,n-l)  matrix.  The 
following  operations  would  occur:  rewind  file,  skip 

TJ_n  9 

forward  — records,  read  m  the  next  record  which  con¬ 
tains  $(n,n-l),  and  read  and  process  the  remaining  records 
sequentially.  The  last  record  contains  $(1,0). 

b.  Indexed  or  random  access  file.  For  this  file 
structure,  each  record  is  identified  by  a  key.  This  key 
is  assigned  at  the  time  the  record  is  filed,  and  it  can  be 
simply  the  integer  of  the  highest  index  in  the  record  (e.g. 
key  =  j  for  the  record  containing  $(kj,kj-l)  ).  As  each 
record  is  filed,  the  record  key  and  the  record  address 
are  stored  in  a  dictionary  for  later  address  retrieval. 
Since  the  key  uniquely  identifies  the  record,  the  file 
chronological  order  does  not  have  to  be  reversed,  just  the 
order  in  which  keys  are  produced  in  file  retrievals.  The 
operations  to  calculate  P(n)  are:  calculate  key  for  the 
record  containing  $(n,n-l),  find  record  address  in  key 
dictionary,  retrieve  record,  process  this  record,  and  con¬ 
tinue  this  cycle  in  a  reverse  chronological  order  until 
the  $(1,0)  matrix  is  used.  The  key  calculation  can  be 
simply  ■  key  =  integer  (-n— ^-~— )  .  For  this  keying  system, 
the  j-th  record  contains  $(kj,kj-l)  through  $(k j-k+1, k j-k) 
data. 
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c.  Indexed-sequential  file.  This  file  structure  is 
designed  for  efficient  operation  in  either  a  random  access 
or  a  sequential  mode.  Each  record  is  assigned  a  key  and 
filed  like  an  indexed  file.  The  records  are  contiguous 
in  storage,  so  once  the  initial  record  is  located  by  key 
search,  the  next  and  subsequent  records  can  be  retrieved  by 
sequential  mode  operations.  The  fact  that  every  record 
retrieved  does  not  require  a  key  search  is  a  significant 
savings  in  some  applications.  To  calculate  P(n) ;  the  file 
containing  $(n,n-l)  is  found  by  calculating  the  key  and 
consulting  the  dictionary;  remaining  records  are  read  into 
core  sequentially  until  the  final  record  containing  0(1,0) 
is  processed.  The  key  for  the  reversed  chronological  file 
ascends  as- time  descends.  This  key  could  be  simply;  key  = 
N+l  -  integer  (--+^~1)  .  The  j-th  record  would  then  contain 
the  0(N- jk+k,N- jk+k-1)  through  0(N- jk+1, N- jk)  data. 

The  sequential  file  structure  was  not  used  because 
the  file  rewind  and  record  skip  operations  were  expected  to 
use  a  prohibitive  amount  of  input-output  time.  Some  simple 
experiments  with  random  access  files  indicated  that  several 
hours  of  input-output  time  would  be  needed  for  one  of  the 
408-point  examples  of  Section  IV.  When  the  indexed- 
sequential  file  structure  was  used  on  these  problems,  the 
input-output  time  required  was  only  1600  seconds  --  less 
than  i  hour.  The  indexed-sequential  file  structure  seems 
uniquely  in  harmony  with  the  data  retrieval  desired  in 
these  analyses.  In  the  studies  performed  during  this 


49 


research,  this  file  structure  afforded  substantial  compu¬ 
tational  costs  savings  over  the  other  structures. 

The  savings  associated  with  indexed-sequential  file 
structure  are  offset  somewhat  by  increased  core  require¬ 
ments  to  accommodate  the  system  software  which  manages 
the  record  storage  and  retrieval.  If  a  sequential  file 
structure  were  used,  these  same  operations  would  be  required 
in  the  program  explicitly,  however.  Additional  savings 
with  indexed-sequential  files  are  gained  by  providing  a 
buffer  in  core  for  several  records  to  come  in  simultaneously. 
For  the  examples  analyzed  in  this  study,  this  increased 
core  storage  requirement  was  judged  to  be  acceptable  in 
order  that  the  input-output  time  requirement  be  acceptably 
low. 

In  summary,  four  distinct  methods  were  used  to  solve 
the  <£>-data  logistics  problem: 

a.  Only  $(i,i-l)  matrices  are  produced  prior  to  the 
P(t)  calculations, 

b.  ^-matrices  are  stored  in  reverse  chronological 
order  for  subsequent  retrieval, 

c.  Several  ^matrices  are  packed  into  each  record  of 
the  data  file  to  decrease  the  number  of  record  retrievals, 
and 

d.  Indexed-sequential  files  are  used  for  efficient 
retrieval. 
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These  methods  essentially  trade-off  computer  core  increases 
for  lower  input- output  time. 

Attention  has  been  focused  on  the  need  to  supply 
$-data  to  the  finite  summations,  but  a  similar  need  for 
r-data  exists  in  the  Q-matrix  evaluations.  The  integrand 
evaluations  can  be  organized  such  that  these  two  data  needs 
are  parallel.  While  the  <$(n,n-l)  to  $>(n,0)  matrices  are 
being  calculated  and  used,  the  pairs  (r  ,r  through 
(rn,rQ)  are  needed  for  the  associated  Q  evaluations.  The 
r-data  can  be  stored  with  the  $-data  by  grouping  r(i)  with 
$(i,i-l).  This  pairing  of  r  and  Ois  convenient  since  the 
coupled  trajectory  and  $  computations  will  naturally  yield 
this  set  at  the  same  time.  Time  data  may  similarly  be 
required  for  At^  calculation.  For  this  study,  the  record 
structure  wass 

(f  -  m  +  1)  , 

r(mk) ,  r(mk-l) ,  ....  r(mk-k+l)  , 

$ (mk,mk-l) ,  $  (mk-l,mk-2) ,  ...,  $ (mk-k+l,mk-k) , 

^mk  ’  \ik- 1 . \ik  -  k+ 1 

N 

where  m  =  1,2,  . . . ,  ^  .  The  record-unique  key  is  the 
N 

integer  ^  -  m  +  1.  In  the  examples  of  Section  IV,  a  four- 
element  representation  of  r  was  used  and  four  (k)  data  sets 
per  record.  This  record  format  requires  345  words  per 
record  versus  the  8?  words  per  record  with  k  *  1.  An 
additional  1300-word  buffer  was  provided  for  the  record 
manager  system  software.  These  computer  core  costs  were 
judged  to  be  modest  in  comparison  to  the  input-output  time 
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savings . 

These  computational  considerations  seem  mundane 
compared  to  the  theoretical  aspects  of  this  work.  But  for 
full-scale  and  realistic  problems,  attention  to  these 
details  will  probably  be  the  deciding  factor  in  whether 
or  not  the  analysis  will  be  done.  The  procedures  outlined 
above  solved  this  problem  within  the  budgetary  constraints 
of  this  research  project.  The  analyses  performed  for 
Sections  IV  through  IX  would  have  been  unacceptably  expen¬ 
sive  had  the  data  logistics  not  been  organized  efficiently. 
These  methods  should  be  useful  in  other  covariance  integral 
approximations.  Whether  to  use  any  or  all  of  these  data 
logistics  methods  should  be  determined  on  a  case-by-case 
basis.  Diverse  factors  which  are  problem  and  computer 
facility  unique  will  certainly  need  to  be  considered. 

The  algorithmic  forms  will  be  developed  in  the  next 
subsections.  The  issues  discussed  here  which  influence  the 
choices  in  algorithmic  form  are: 

a.  Decision  to  produce  P(n)  recursively,  and 

b.  The  efficiencies  of  reverse  chronological  sum¬ 
mation. 

These  factors  should  be  recalled  in  the  algorithmic  develop¬ 
ments  which  follow. 

3.2  Two  Dimensional  Direct  Integral  Approximations 

With  these  data  handling  issues  concluded,  the  direct 
approximation  of  (  2$  can  be  addressed.  The  task  is  to 
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produce  the  entire  sequence  of  the  time  history  [P(  0), 

P(  1 ) ,  ....  P(  K)},  not  just  one  element.  So,  a  sequential 
algorithm  is  suggested. 

To  begin  this  development,  a  general  discrete  approxi¬ 
mation  of  Eq  (2?)  can  be  written  as 

P(  n)  =  2  2  $(n, i)G( i) Q( i, j)GT( j)$T(n, j)  • 

i=0  j=0 

•  Atj  Ati  sn(i,j)  (30) 

where 

1)  At.  and  At.  are  time  intervals  associated  with  the 

J 

i-th  and  j-th  data  points. 

2)  sn(i»j)  is  a  discrete  weighting  function  which 
depends  on  the  index  n.  A  different  function  s  (  •  ,  • )  will 
be  defined  for  each  integration  rule. 

Note  that  for  proper  weighting 

2  2  s  (i,  j)  At.  At.  =  (t  -t  )2  (33) 

i=0  j=0  n  1  J  n  0 

The  simplest  choices  for  sn  functions  correspond  to 
rectangular  and  trapezoidal  integration  rules.  Others 
can  certainly  be  developed,  but  these  simple  rules  will 
yield  the  least  cumbersome  recursion  relationships.  The 
improvements  from  re c tangular^'t'o  trapezoidal  will  give 
insight  into  the  incentive  for  higher  order  development. 

3.2.1  Rectangular  Integration  Approximation  Algorithm 
The  rectangular  integration  approximation  gives 
equal  weighting  to  all  evaluation  points.  The  formulation 
here  is  a  slight  modification  of  the  usual  rectangular 
rule,  but  the  essential  character  of  the  method  is  retained. 
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In  one  dimensional  integration  approximation,  the  integral 
is  calculated  using  n  integration  steps  and  n  integrand 
evaluations.  In  the  following  development  n+1  integrand 
evaluations  are  assumed  to  bracket  the  integration  inter¬ 
val.  Each  integrand  evaluation  in  the  finite  sum  is 
adjusted  by  a  weighting  factor  of  n/n+1  to  account  for 
this  change.  In  this  subsection,  this  simple  one-dimensional 
integral  approximation  rule  will  be  extended  to  two  dimen¬ 
sional  (double)  integrals.  This  rule  will  be  interpreted 
with  (30)  to  form  a  recursion  which  yields  the  sequence 

£P(  i  ) 3- 

The  one-dimensional  weighting  of  integrand  evalua¬ 
tions  can  be  expressed  as  a  sequence  of  functions.  The 
generic  element  is  defined  for  n  >  1  by 


6nU>  = 


i  <  0  or  i  >  n 
0  <  i  <  n 


(34) 


This  sequence  weighting  function  gives  the  weight  assigned 
to  the  i-th  integrand  evaluation  for  i  =  0,1,2,  ....  n. 

The  extension  to  two  dimensions  is  simply 

sn(ili)  '  6n(i)  6n(j> 


for  both  0<i<n  and  0<j<n 
otherwise 


(35) 


n  n  ? 

Note  that  _2  2  s  (i,j)  =  n  ,  so  this  constant  weighting 

1=0  i=0 
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A  con- 


n 

will  satisfy  (33)  as  long  as  2  At.  =  t  -  t  . 

i=0  i  n  o 

stant  step  size  was  used  in  this  work  with  At^  =  At  = 

(tn-tQ)/n  for  0  <  i  <  n.  A  more  general  variable  step 

size  approach  might  prove  more  efficient  on  missions  with 

a  wide  range  of  dynamics.  This  generalization  is  not 

required  for  the  verifications  and  demonstrations  of  this 

work,  but  a  suggestion  for  future  research  along  this  line 

is  given  in  Section  XI. 

With  s  (i, j)  defined  by  (35),  the  direct  approximation 
of  (27)  could  proceed  according  to  (30)  for  each  value  of 
n.  A  recursion  is  more  efficient  when  a  sequence  of  answers 
is  desired.  The  solution  of  the  homogeneous  form  of  (27) 
is  well-known  [Ref  15  ••  Equation  (4-114)  on  page  165  with 
Q(r)  =  0].  With  some  algebraic  manipulations,  the  term, 
$(n,n-l)  P(n-l)  $  (n,n-l) ,  appears  within  the  P(n)  summa¬ 
tion  and  forms  the  basis  for  the  recursive  relationship. 

To  develop  this  recursive  relationship,  write  (30) 


where 


P(n)  =  +  S2  +  (36) 

a 

S,  *  2  2  $(n,i)G(i)Q(i,j)GT(j)$T(n,j)  At.  At.s  (i, 3) 

x  i=0  3=0  1  j  n 


J  \  riT  t  .  \  *T  , 


>0  =  2  $(n,n)G(n)Q(n,  o)G  (j)$  (n,j)  At  At.s  (n, 3) 

c  3=0  n  3  n 

(38) 
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S,  =  z  $(n,i)G(i)Q(i,n)G  (n)<2>  (n,n)  At.  At  s  (i,n) 

3  i=0  xnn 

(39) 

=  $>(n,n)G(n)Q(n,n)GT(n)$T(n,n)  At^sn(n,n)  (40) 
Some  obvious  simplifications  can  be  made  since  $(n,n'  =  I. 
A  decomposition  of  (37)  should  yield  P(n-l),  and  the  fac¬ 
tored  form  is  obtained  using 

$(n,i)  =  $(n,n-l)<£>(n-l, i)  (41) 

which  yields 


S,  =  <$(n,n-l)(  2  2  $(n-l,i)G(i)Q(i,  j)GT(  j)$T(n-l,  j) 

1  l  i=0  j=0 

Ati  At^sn(i,j)l  •$T(n,n-l) 

J  (42) 


Note  that  the  indices  of  summation  in  (42)  are  such  that 

2  2 

sn  i^’j)  is  (n-1)  /n  throughout  the  ranges.  Since 
2  2 

sn(i,j)  is  n  /(n+1)  ,  a  replacement  can  be  made  using 


sn(l,j)  =  (n2^1)2  sn-l(l’j) 


(43) 


for  both  0  <  i  <  n-1  and  0  <  j  <  n-1.  With  this  relation¬ 
ship 


s  = - g - ?  $(n,n-l)(  2  2  $(n-l,  i)G(  i) Q(  i,  j)GT(  j) 

1  (n2-ir  l  i=0  j=0 


$T(n-l,j)  At^  At^sn_1(i, j) |  $T(n,n-l) 

-  - V—o  $(n, n-1)  P(n-l)$T(nfn-l)  (44) 

(n2-l)2 

The  homogeneous  solution  is,  in  this  case,  modified  by  a 
factor  which  reflects  the  change  in  weighting  from  tn  ^  to 
the  tn  evaluation  point.  S2,  S^»  and  represent  the 
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change  in  P  over  the  t  to  tn  interval  due  to  new  dis¬ 
turbances  entering  the  problem.  In  that  sense,  S2»  S^, 
and  represent  the  current  time  effects  while  S-^  repre¬ 
sents  the  past  time  effects.  This  type  of  interpretation 
is  important  later  in  applying  these  algorithms  to  more 
complex  trajectory  and  correlation  models  than  those  which 
follow  in  this  section. 

Equation  (44)  establishes  the  basic  recursive  rela¬ 
tionship  with  (3 6).  Some  additional  interpretations  will 
ease  the  computational  burdens  associated  with  S2  and  S^. 
First,  from  the  definition  of  Q,  (28),  it  should  be  clear 
that 

Q(i,  j)  =  QT(j,i)  W) 

Using  this  relationship,  it  follows  that 

S3  =  (46) 

This  relationship  affords  a  considerable  computational 
savings  since  these  summations  are  the  major  part  of  the 
calculations  required  on  each  recursion. 

These  "current  time"  contributors  to  P(n)  have  a  strong 
similarity  to  the  D  function  of  Nested  Integrals,  (29). 

At  this  point,  the  S2  and  terms  will  be  redefined  in  a 
way  which  will  enhance  one's  appreciation  of  the  likeness. 
Define 


n-1 


Now,  equation  (41)  can  be  restated  as 

„4  m 

P(n)  =  — = - P  $(n,n-l)P(n-l)$(n,n-l)  +  D  (n)  +  D*(n) 

(n2-l) 2  a  a 

+  Db(n)  (49} 

The  recursion  is  defined  by  (4?),  (48),  and  (49),  but  the 
summation  of  (4?)  must  be  properly  interpreted  to  yield 
the  reverse  chronological  form.  To  make  this  point  expli- 


D  (n)  =  G(n)  2  Q(n,n-i)GT(n-i)$T(n,n-i)  At  At  .s(n,n-i) 

<*  •  n  n  11“  X  II 

(47a) 

Slimming  in  the  order  prescribed  by  (47a)  requires  $(n,n-l), 
$(n,n-2) ,  ....  $(n,l) ,  and  $(n,0)  in  that  order.  This 
arrangement  allows  the  data  logistics  economies  discussed 
in  subsection  3.1.  Figure  8  presents  the  algorithmic 
structure  which  executes  the  summation  of  (47a).  Note  that 
at  the  top  of  the  DO  loop  the  local  variable  $Temp  =  $(n,j), 
and  note  that  the  multiplication  at  the  bottom  of  the  loop 
requires  only  $(j,j-l)  as  new  data.  Since  there  is  no 
$(0,-1),  the  rQ  evaluation  point  is  treated  separately.  The 
READ  operations  should  not  be  taken  literally.  Recall 
that  several  data  sets  may  be  packed  in  one  record;  therefore, 
READ  means  "locate"  which  may  or  may  not  require  an  actual 
READ  operation.  The  sequential  data  retrieval  within  the 
loop  is  initialized  by  the  r(n)  data  retrieve  upon  entry 
to  the  algorithm.  Again,  these  data  handling  techniques 
are  not  the  only  solution,  but  algorithms  formed  in  this 
manner  allowed  this  work  to  proceed. 
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The  recursion  of  (49)  applies  for  n  >  2  because  of  the 
obvious  singularity  at  n=l.  Since  P(0)  is  identically 
zero  by  previous  assumptions,  one  might  make  the  mistake 
of  using  the  last  three  terms  of  (49)  to  produce  P(l). 

This  would  neglect  the  Db(0)  term  which  occurs  in  (30) 
evaluated  for  n=l.  This  term  adds  the  effects  of  disturb¬ 
ances  near  the  start  of  the  mission.  The  first  step,  then, 
can  be  calculated  by  executing  (30)  for  n=l. 

P(l)  =  <£(l,0)Db(0)$T(l,0)  +  Da(l)  +  d£(1)  +  Db(l) 

(50) 

The  algorithm  is  completely  defined  by  (47a),  (48), 
and  (49)  for  the  recursion  and  by  (47a),  (48),  and  (50) 
for  the  first  step.  The  algorithmic  procedure  is  straight¬ 
forward  except  for  the  D_(n)  computation  shown  in  Figure  8. 

CL 

For  each  iteration  after  t^,  Da(n)  and  Db(n)  are  computed. 
Then,  P(n)  is  computed  by  (49)  using  the  $>(n,n~l)  which 
is  on  file. 

3.2.2  Trapezoidal  Integration  Approximation  Algorithm 

The  Trapezoidal  algorithm  will  next  be  developed  in  a 
manner  similar  to  the  Rectangular  development.  First,  a 
weighting  function  will  be  specified,  and  then,  recursion 
and  initialization  relationships  will  be  formed. 

The  weighting  function  for  one-dimensional  integra¬ 
tion  approximation  by  the  trapezoidal  rule  is 
JO  i<0  or  i>n 

6n(i)  =  <  i  i  =  0  or  i  =  n  (51) 

1  1  <  i  <  n-1 
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Extending  this  to  two  dimensions,  define 


sn(i,j)  = 


sn(i, j) 


6n(1)6n<j) 

(i,j)  9  ( i<0  or  j<0  or  i>n  or  i>n) 

(i»  j)  «  C(0.0) ,(0,n) ,(n,0) , (n,n)3 

( i . 3 )  «  [(O.k) , (k,0) , (n,k) , (k,n)  9  l<k<n-l] 
(i,j)  ?  (1  <  i  <  n-1  and  1  <  j  <  n-1) 

(52) 


The  rule  expressed  by  (51)  applies  to  equal  subintervals 
[Ref  16 s 144]  which  will  apply  in  the  analyses  following. 
This  form  can  be  used  for  non-uniform  step  sizes  but  care 
must  be  taken  to  satisfy  the  constraint  of  (33) . 

An  alternative  for  nonuniform  step  sizes  is  to  give 
unity  weighting  to  all  integrand  evaluations  and  to  com¬ 
pute  step  size  as  a  function  of  n  by 


i 


ri(ti  -  to) 

i 

=  0 

ti(n)  =< 

1  ^('ti+l"ti-l) 

1 

<  i  <  n-1 

(53) 

^(tn  -  tn-l) 

i 

=  n 

For  this  development  sn(i,j)  is  retained  since  it  aids  the 
development  of  a  recursive  relationship. 

In  the  development  of  the  Rectangular  recursion,  use 
was  made  of  the  constant  ratio  between  sn(i,j)  and 
sn_1(i,j)  over  a  restricted  range  of  (i,j).  In  the 
Trapezoidal  case,  a  similar  useful  property  exists  only 
this  time  it  is  the  difference  of  the  two  weighting  func¬ 
tions.  Note  that  sn(i,j)  -  sn_p(i»j)  Is  zero  except  for 
the  cases  of  either  i  or  j  equal  to  either  n  or  n-1. 


sn(i,j)  -  sn_i(i>0) 

The  weighting  functions  are  the  same  over  indices  ranging 
from  0  to  n-2.  Over  the  range  of  0  to  n-2  of  one  index 
with  the  other  fixed  at  n-1,  sn_i(i»0)  is  exactly  one-half 
of  sn(i,j).  The  n-1  corner  element  shows  a  4:1  ratio  of 
sn(n-l,n-l) : sn  1(n-l,n-l) .  With  these  relationships  in 
mind,  the  summation  of  (30)  is  partitioned  by 

P(n)  =  Sx  +  S2  +  S3  +  S4  +  S5  +  S6  +  S?  (54) 

where  for  this  Trapezoidal  development  the  0  to  n-2  sub¬ 
block  gives 

S,  =  2  2  $(n, i)G(i)Q(i, j)GT( j)$T(n, j)  At.  At.s  (i,j) 

1  i=0  j=0  1  J  n 

(55) 

The  elements  of  the  sum  with  n-1  index  in  one  position  give 
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n  2 

S9  =  2  $(n,n-l)G(n-l)Q(n-l,i)GT(i)$T(n,i)  At,  At  , 

^  i=0 

sn(n-l,i)  (56) 

Using  (45)  again  the  next  term  can  be  stated  as 

S3  =  Sg  \  (57) 

And  the  (n-l,n-l)  corner  element  is  given  by 

S4  =  $(n,n-l)G(n-l)Q(n-l,n-l)GT(n-l)2>T(n,n-l)  At^_1 

sn(n-l,n-l)  (58) 

The  elements  with  n  as  the  first  index  are  included  in 
n  X 

s,  =  2  3>(n,n)G(n)Q(n,i)G^(i)<!>T(n,i)  At,  At  s  (n,i)  (59) 

i=0  inn 

Using  ( 45)  ,  the  other  n-index  elements  are  given  by 


The  (n,n)  corner  element  is 

S?  =  $(n1n)G(n)Q(n,n)GT(n)$'r(nf  n)  At^sn(n,n)  (61) 

The  sn(n-l,i)  of  (56)  can  be  replaced  by  2sn_1(n-l,i) 
using  the  relationship  discussed  previously.  Similarly, 
sn(n-l,n-l)  in  (58)  can  be  replaced  by  4sn_2(n-l,n-l) . 
Using  these  changes  and  factoring  out  $(n,n-l) ,  yields 

S1  +  |S2  +  |S2  +  iS  4  =  $(n,n-l)P(n-l)$T(n,n-l)  (62) 

With  this  result,  (5*0  can  be  restated  as 
P(n)  =  £>(n,n-l)P(n-l)$T(n,n-l)  +  |S2  +  £s£  +  OA)S4 

+  s5  +  +  s7  (63) 

With  the  definitions, 

XI  1 

D  (n)  =  S,  =  2  G(n)Q(n,  i)GT(i)$T(n,  i)  At,  At  s  (i,j) 
a  5  i=0  1  n  n  (64) 
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(65) 


Db(n)  =  Sy  =  G(n)Q(n,n)GT(n)  At^sn(n,n) 

the  following  identifications  can  be  made 

iS2  =  $(n,n-l)Da(n-l)$T(n,n-l)  (66) 

(3A)S4  =  3^(n,n-l)Db(n-l)$T(n,n-l)  (6?) 

Drawing  all  of  these  terms  together  yields  the  desired 
recursion 

P(n)  =  $(n,n-l)  {P(n-l)  +  Da(n-1)  +  D^(n-l)  +  3Db(n-l)} 

$T(n,n-l)  +  Da(n)  +  Da(n)  +  Db(n) 

(  68) 

Again  the  main  computational  burden  is  the  D0  compu- 

a 

tation  and  the  reverse  chronological  order  is  implemented 

by 

Da(n)  =  i  G(n)Q(n,n-i)GT(n-i)$T(n,n-i)  Atfi  Atn_isn(n»n~i) 

(64a) 

The  Trapezoidal  algorithm  recursion  is  defined  by  (64a), 
(65)1  and  (68).  The  P(l)  calculation  of  the  Rectangular 
method  is  used  for  the  first  step  since  for  n=l  the  two 
methods  are  equivalent.  Figure  8  is  a  valid  description 
of  the  Trapezoidal  algorithm  except  the  past  value  of  D 

a. 

must  be  saved. 

3*  3  Nested  Integrals  Approximation  Algorithm 

As  an  alternative  to  these  two  direct  approximations 
of  (  27) ,  an  algorithm  based  on  the  sequential  solution  of 
two  one-dimensional  integrals  is  suggested.  This  nested 
pair  of  integrals  results  when  one  takes  the  derivative  of 


* 


V 


(  2$  with  respect  to  time,  t.  The  resulting  algorithm  is  . 

similar  to  those  in  the  direct  integral  approximations  by  * 

rectangular  and  trapezoidal  rules.  This  new  form  permits 
the  use  of  powerful  ordinary  differential  equation  tech- 
niques  to  produce  the  final  result.  • 

The  first  step  in  the  development  is  an  application 
of  Liebnitz'  rule  to  (27): 

P(t)  =  G(t)  Q(t,q)  GT(q)$T(t,q)dq 

o 

t  _  m 

+  /+  $(t, p)G(p)Q[r(p)  ,r(t)  ]dp  Gx(t)  (69) 

xo 

+  F(t)  P(t)  +  P(t)FT(t) 

Defining 

D(t)  =  f  Q(tfp)GT(p)^T(t,p)dp  (29b) 

0 

and  noting  from  (45)  that 

Q(t,p)  =  QT(p,t)  (**5a) 

results  in 

P(t)  =  G(t)D(t)  +  DT(t)GT(t)  +  P(t)FT(t)  +  F(t)P(t)  U9a) 

The  differential  equation  in  (29a)  together  with  its  associ¬ 
ated  equation  (29b)  comprise  the  nested  integrals  pair. 

To  solve  these  equations,  first,  approximate  a  solu¬ 
tion  to  (29b)  for  D(t)  over  the  full  span  of  mission  time 
for  which  P(t)  is  required.  A  completely  new  integral  j 

occurs  for  each  t  since  r(tn)  occurs  as  an  argument  to 
Q(  •  ,  • )  for  the  first  time  in  the  n-th  integral  of  the  1 

sequence.  j 
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D(  n)  =  sln  QCr(tn),r(p)]GT(p)$!r(tn,p)dp  (70) 

o 

Since  [D(  1),D(  2) . D(  N  )  }  must  be  produced,  let  the 

approximation  take  the  form  of  a  simple  quadrature  formula. 
Using  the  reverse  chronological  form, 

D(n)  =  2  Q[r(  n)  ,r(  n-i  )]GT(  n-i  )$T(  n  ,  n-i  ) 

i=0 


*  Atn-isn(n"i) 


(71) 


One  could  develop  sn(i)  for  an  involved  quadrature  form. 

To  keep  the  computations  simple,  lower  the  time-step  incre¬ 
ments  to  reduce  integral  approximation  error  to  an  accept¬ 
able  level.  So,  define 

r 


sn(i)  - 


1 

2 


i 


i  =  0 

1  <  i  <  n-1 
i  =  n 


►(72) 


corresponding  to  trapezoidal  integration.  Then  evaluating 
(71)  for  n=l, 2,  ...,  N  will  yield  a  sequence  [D(  0) ,D(  1), 
...,  D(  N  ) }  which  can  be  used  in  solving  (29a). 

One  can  treat  (29a)  with  a  multitude  of  ordinary 
differential  equation  methods.  For  the  following  numeri¬ 
cal  examples,  a  predictor-corrector  £Ref  14J  was  employed, 


*  As  previously  discussed,  this  weighting  could  also  be 
viewed  as  an  alteration  or  redefinition  of  At^.  The  form 

presented  here  makes  such  weighting  explicit  for  clarity. 
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*  because  P  is  a  function  of  P.  Since  evaluations  of  (29b) 
are  required  at  points  other  than  CtQ,t^,  ....  t^},  a 
simple  linear  interpolator  was  used  on  D(t) ,  F(t),  and 
G(t).  The  overall  Nested  Integrals  algorithm  is  a  hybrid 
integration  method  using  trapezoidal  quadrature  for  D(t) 
and  predictor-corrector  methods  for  P(t). 

Here,  one  should  note  that  the  Nested  Integrals  algo¬ 
rithm  requires  F( t)  data  explicitly,  whereas  the  direct 
integral  approximations  did  not.  The  data  flow  is  some¬ 
what  altered  by  this  need.  The  flow,  including  the  indexed 
sequential  file  step  (subsection  3.1),  is  diagrammed  in 
Figure  9.  The  D(n)  computation  is  quite  similar  to  the 
D„(n)  computation  portrayed  in  Figure  8.  Obvious  changes 
are  the  (n,n)  point  and  the  deletion  of  G(n)  and  Atn  terms. 

3.4  Comparison  Rationale 

One  can  speculate  about  the  relative  merits  of  these 
three  algorithms;  Rectangular,  Trapezoidal,  and  Nested 
Integrals.  The  spatial  spectrum  of  the  anomalous  gravity 
correlation  functions  and  the  temporal  spectrum  of  the 
inertial  navigation  error  propagation  are  the  principal 
concerns  in  selecting  step  size  for  any  method.  Nevertheless 
all  three  methods  are  expected  [Ref  17:382-384]  to  yield 


*  The  specific  predictor-corrector  was  variable-order, 
variable-step-size  Adams-Pecce  algorithm  described  in  Ref 
14.  Other  predictor-corrector  formulations  of  comparable 
order  should  give  comparable  results. 
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Figure  9.  Nested  Integrals  Information  Flow 


acceptable  accuracy  when  both  the  products  of  t^  with 
F(t)  eigenvalues  are  small  compared  to  1  and  |r(t^)  - 
r(t^  2.)  I  is  small  compared  to  correlation  distances  in 
Q(.,.).  The  Trapezoidal  algorithm  should  be  somewhat  more 
accurate  than  the  Rectangular,  because  a  higher  order 
interpolating  polynomial  represents  the  function  between 
evaluation  points. *  The  amount  of  improvement  here  will 
indicate  the  productivity  of  more  complex  direct  approxi¬ 
mations.  The  predictor-corrector  P(t)  computation  in  the 
Nested  Integrals  algorithm  should  give  an  accuracy  advan¬ 
tage  but  will  undoubtedly  cost  in  computer  storage  and 
execution  time. 

Accuracy  and  efficiency  will  be  used  to  select  one 
of  the  alternatives  for  further  development  and  study. 

To  quantify  the  differences  for  this  selection,  a  simple 
test  case  is  proposed.  This  example  should  be  in  the  form 
of  the  problem  type  that  these  algorithms  are  intended  to 
solve,  but  it  should  be  simple  enough  for  an  indisputable 
comparison  to  the  true  solution. 

3.5  Undamped  Schuler  Loop  Driven  by  Exponentially  Correlated 

Noise 

A  simple  model  of  an  inertial  system  horizontal  channel 
is  the  Schuler  loop.  In  its  purest  form,  this  model  is  an 


*  Hildebrand  £  Ref  18  *  93—  9^1  has  shown  that  for  interpo¬ 
lating  functions  of  order  2n  and  2n+l  the  approximation 
error  is  bounded  by  a  term  involving  A  t^n+ 3, 
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undamped,  second-order  feedback  loop.  The  feedback  signal 
represents,  to  first  order,  the  corrective  influence  of 
the  gravity  model,  given  position  errors  [Ref  2  «  Appendix 
B],  The  feedback  weighting  is  the  squared  Schuler  fre¬ 
quency,  (l)g  =  f“  . 


6g(t) 


As  diagrammed  above,  this  loop  is  driven  by  the  gravity 

model  error  term  6g  which  corresponds  to  the  vertical 

deflection  in  the  +x  direction  multiplied  by  g. 

This  double  integration  with  negative  feedback  forms 

a  pure  oscillator  with  fundamental  frequency  (j  . 

s 

Defining  state  variables  as  integrator  outputs  yields 


f  6*”)  fO  ll  f  6x^  fo\ 

-0)g2  0  \6v/+\l/ 


6g 


(73) 


For  the  algorithms,  identify 


From  ( 74) ,  one  can  derive 
$>(t,p)=$>(t-p)  = 


coso)s(t-p) 


-Ws  sind)s(t-p) 


—  sinu8(t-p) 
s 

COS(i)s(t-p) 


(78) 


Levine  and  Gelb  [Ref  ll]  suggest  a  simple  exponential 
correlation  for  such  a  6g.  The  correlation  is  spatial 


since  the  gravity  model  error  is  a  spatial  function. 

|x„-x,| 


(g)^6g(x1)  6g(x2)J  =  a2« 


(79) 


2  - 


where  a  is  the  variance  and  d  is  the  correlation  distance. 

From  (79) ,  form  the  required  gravity  correlation  func¬ 
tion 


IxttgJ-xCt-^  t 


QLlltt-^  ,r(t2)  ]  =  Q[x(  tx)  ,x(t2)  ]  =  o2  e 


(80) 


Only  the  design  mission  remains  to  be  specified. 

Again,  simplicity  is  desired,  so  a  constant  velocity  trajec¬ 
tory  is  selected 

x(t)  =  V  •  t  (81) 

Now, 

|x(  t2)-x(  tx)  |  =  V  •  Itg-t-jJ  (82) 

in  (80)  above.  For  the  numerical  example,  data  is  specified 
by  Table  II.  These  data  correspond  to  values  for  horizontal 
gravity  modeling  errors  [Ref  12]  and  a  bomber  cruise  mis¬ 
sion.  The  intent  is  to  form  only  a  representative  case, 
so  precision  in  these  data  is  not  an  issue. 
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Table  II 


Schuler  Cases  Data 

Quantity _ .  Symbol  Value,  Units 


Initial  time 
Velocity 

Disturbance  Variance* 
Correlation  Distance 
Schuler  Frequency 


t  0  sec 

o 

V  615  ft/sec 

a2  1273  (mgal)2 

d  20  n  m 

0_  1.24xl0~3 

s  sec 


*  Corresponds  to  a  36  pg  root-mean-square  which 
is  equivalent  to  7-5  arc  sec  rms. 

1  mgal  =  10-3  cm/sec2  =  10"3  m/sec2. 


3.5.I  Closed  Form  Solution 

The  information  to  perform  the  integral  approxima¬ 
tions  is  now  completely  specified.  Now  P(t)  must  be  pro¬ 
duced  in  closed  form  to  study  and  compare  algorithm  errors. 

The  problem  can  be  solved  using  a  straightforward 
approach  based  on  stochastic  linear  system  theory.  The 
model  for  gravity  error,  6g,  is  the  output  of  a  linear 
filter  driven  by  white,  gaussian  noise.  The  block  diagram 


where 

<^Eqx(xl)qx(x2^  =  ^d"  6(x2“X1)  *83) 

Also,  the  filter  is  modeled  at  steady  state, 

<f[6gx(x)2]  =  c2  (84) 

This  spatial  model  is  converted  to  temporal  by  the  velocity. 
Simply  replace  dx  by  V*dt,  and  the  block  diagram  can  be 
converted  to 


q(t) 


Jdt 


6g(t) 


3 


with 


Q  =  X 

p  d 


(85) 


The  noise  process  satisfies 

(fCqfVqCtg)]  =  2 3ct2  6(t2-tx)  (85) 

Combining  the  Schuler  loop  with  this  model  yields  one 
linear  system  driven  by  q(t).  An  augmented  state  equation 
can  now  be  formed 


6x 

6v 

6g 


0 

2 

-ws 

1 

0 

1 

1 

1 

1 

1 

1 . . . 

1 

O  iH  1 
1 

0 

0 

t 

1 

1 

-3  J 

6x 

6v 

6g 


0 

0 


q(t) 


x„(t) ,  F„ ,  and  GQ  are  given  by  ( 7  )  through  ( 10)  . 

"**a  a  cl 

Pa(t)  =  F_P_(t)  +  P_(t)  F^  +  G  Q  G 


a  g  a 


where  Qg  comes  from  (  6 )  as 

«g  -  [23  O2] 


(87) 

Then, 

(14) 

(88) 
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Define 


'Pll 

P12 

P13 

Pa  = 

P12 

P22 

P 

23 

lP13 

p 

23 

p 

33 

note  that 

'0 

0 

!  o  ' 

i 

pa<°> 

- 

0 

0 

!  o 

_0 

0 

■f - 

1  °2. 

Now  define 

h  =  0/(jJs 

1 


A  = 


02+Ws2 


Then  solving  (14)  yields 


(89) 


(90) 


(91) 

(92) 


pu<t> 

=  G-^Jht 

ws 

~  -7T—  [h  sin 2(o 
2u)g 

t+1  -  cos2a)  t]  + 
s  s 

2A(I)g[l-e"pt(h 

sinwot  +  cos(0_t)J} 

S  o’ 

(93) 

r12<*> 

2 

=  Q^p{|[h(l-cos2(Jost)-sin2wst]  +  e-^"1  sinu)gt} 
s 

(94) 

P22<*> 

=  a2A[8t 

+  |(hsina)  t  - 

S 

cos2(i)st+l)  - 

2(i)  RACh+ 

D 

e_^^(  sin(ost  -  h  cosa)s 

t)]} 

(95  ) 

V'l 

=  a2A£l-e 

’^^[cosai  t  +  h 
s 

sind)gt]} 

(96  ) 

p23(t) 

=  a2A[3+e 

“^[-3  COSGJgt 

+  ws  sinu)gt]} 

(97  ) 

P33(t) 

=  a2 

(98  ) 
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Now  (93) »  (94),  and  (95  )  give  the  true  answers  against 
which  the  algorithms’  performance  can  be  judged.  Note 
that  the  position  and  velocity  variance  terms  have  ramp 
components.  The  results  (  96)  and  (97  ),  incidentally, 
give  the  correlation  of  position  and  velocity  errors  with 
present  gravity  error.  Equation  (  98  )  merely  confirms  the 
problem  specification  in  (84). 

3.5.2  Numerical  Comparisons 

Trapezoidal  and  Rectangular  algorithm  solutions  were 
formed  using  (75) .  (78),  (80),  and  (81).  The  Nested 
Integrals  algorithm  also  required  (74) .  For  each  of  the 
three  methods,  a  fixed  integration  step  size  was  used 
throughout  an  analysis  case.  Six  different  integration 
step  size  cases  were  considered:  3*?5»  7-5.  15.  30,  60, 
and  120  seconds.  The  errors  in  P(t)  elements  for  each 
of  these  18  cases  are  plotted  against  a  common  scale  for 
comparisons.  The  position  and  velocity  variance  results 
are  presented  in  Figures  C-l  through  C-6  .  The  position- 
velocity  covariance  results  are  presented  in  Figures  C-7 
through  C-9.  All  of  these  graphical  results  are  located 
in  Appendix  C.  These  graphs  verify  that  each  approxima¬ 
tion  technique  is  accurate,  given  a  small  enough  integra¬ 
tion  step  size,  At.  The  variance  approximation  errors 
grow  with  time,  but  the  variances  being  approximated  have 
ramp  terms  as  seen  in  (93)  and  (95  )  above.  A  better 
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perspective  is  gained  if  these  errors  are  expressed  in  a 
percent  of  the  true  value.  The  results  for  percent  posi¬ 
tion  variance  error  are  presented  in  Figures  C-10,  C-ll, 
and  C-12  for  the  three  candidate  methods.  The  only  signi¬ 
ficant  errors,  from  a  practical  viewpoint,  are  the  initial¬ 
ization  errors  associated  with  the  Rectangular  and  Trapezoidal 
algorithms.  The  absolute  error  may  grow,  but  the  percent 
errors  fall  to  clearly  acceptable  levels  for  small  ht. 

Obviously,  for  all  methods,  some  catastrophic  failure 
awaits  the  unwary  user  who  increases  At  beyond  the  largest 
value  (120  seconds)  used  in  this  study.  The  explanation 
is  found  in  applying  Shannon’s  sampling  theorem  to  the 
choice  of  integration  step  size,  and  a  digression  on  this 
point  is  needed  before  the  accuracy  comparison  is  con¬ 
ducted. 

Shannon's  sampling  theorem  dictates  that  be  at 
least  twice  the  highest  frequency  [Ref  15:2953  which 
affects  the  integrand  of  (27).  This  integrand  is  affected 
by  the  error  propagation  model,  the  trajectory,  and  the 
correlation  function.  The  trajectory  is  not  reflected 
in  the  simple  error  propagation  model,  but  it  affects  the 
correlation  function  through  its  arguments.  The  frequency 
characteristics  of  the  propagation  model  and  of  the  trajec¬ 
tory-driven  error  of  propagation  function  should  therefore 
be  investigated. 

The  Schuler  loop  of  this  example  acts  as  a  low  pass 
filter  to  the  gravity  disturbance  input.  The  filter 


dynamics  are  adequately  modeled  when  At  is  less  than  one- 
tenth  the  Schuler  period.  This  imposes  the  constraint  that 
integration  step  size  be  less  than  8.4  minutes  -  clearly 
met  by  all  values  of  At  considered.  The  input  correlation 
function  then,  must,  be  the  key  to  the  problem  for  the 
greater  At  values. 

The  correlation  function  is  spatial,  but,  as  seen, 
velocity  converts  this  process  into  the  time  domain.  The 
correlation  time, 

tc  =  |  =  dA.  (99) 

is  198  seconds  in  this  case.  Adequate  sampling  dictates 

that  At  be  less  than  99  seconds  (it),  and  common  prac- 

c 

tice  [Ref  15:295]  indicates  one  should  select  At  as  low  as 
20  seconds  (tc/l0).  A  99-second  step  size  corresponds  to 
4.7  or  the  logarithmic  integration  step  size  scale  in 
Figures  C-l  through  C-12,  and  a  20-second  step  size  corre¬ 
sponds  to  2.4.  The  sampling  theorem  violation  for  the 
120-second  (5  on  the  logarithmic  scale),  case  explains  the 
significant  increase  in  error  for  that  case.  The  60 
second  (4  on  the  logarithmic  scale)  case  is  near  this 
violation.  The  results  for  60  seconds  and  120  seconds, 
therefore,  do  not  represent  the  algorithm  performance  and 
should  be  eliminated  prior  to  comparing  the  results  for 
the  three  alternative  algorithms. 
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The  results  for  position  and  velocity  variance  error 
were  replotted,  then,  for  the  four  lowest  integration  step 
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sizes.  These  results  are  shown  in  composite  form  in 
Figures  C-13  and  C-14.  This  format  facilitates  the  direct 
comparison  of  approximation  results. 

Schuler  rate  undulations  stand  out  in  the  Rectangular 
algorithm  results.  The  lack  of  predominant  Schuler  rate 
errors  in  the  Nested  Integrals  and  Trapezoidal  must  be 
attributed  to  the  higher  order  representation  of  the  inte¬ 
grand.  Recall  that  Nested  Integrals  has  a  trapezoidal  rule 
approximation  for  D(t),  so  the  similarity  of  these  results 
is  not  so  surprising.  The  generally  lower  and  less  per¬ 
turbed  Nested  Integrals  and  Trapezoidal  computation  error 
gives  these  methods  a  decided  edge  over  the  Rectangular 
alternative.  The  Trapezoidal  maximum  error  (sup  norm) 
in  velocity  variance  is  slightly  lower  than  Nested  Integrals. 
The  Nested  Integrals  is  slightly  better  in  the  position 
variance  calculation  by  the  same  judgment  criteria. 

Accuracy  is  not  the  only  concern.  The  computational 
burden  should  also  be  considered  in  the  final  selection. 

The  input-output  computer  costs  were  practically  the  same 
for  all  three  methods.  The  computer  memory  required  for 
Nested  Integrals  was  only  2000  words  more  than  the  alter¬ 
natives,  which  is  not  a  major  consideration.  The  computa¬ 
tion  time  did  vary,  and  these  times,  normalized  to  the 
smallest  value,  are  given  in  Table  IH. 


Table  III 


At 

SEC 

Normalized  Computation  Time 

NESTED 

RECTANGULAR  TRAPEZOIDAL  INTEGRALS 

3-75 

1485 

1556 

449 

7-5 

370 

387 

114 

15- 

94 

100 

29 

30. 

24 

25 

8 

6o. 

6.3 

6-5 

2.4 

120. 

1.9 

1-9 

1.0 

Note  that  for  all  three  methods,  halving  the  integra¬ 
tion  step  size  doubles  the  number  of  points,  and  increases 

p 

the  computation  time  by  approximately  2  .  This  reflects 
the  underlying  two-dimenstional  integration  which  is 
being  approximated. 

The  surprising  result  of  the  computation  time  compari¬ 
son  is  that  the  Nested  Integral  method  is  more  efficient 
than  the  other  methods  by  a  ratio  of  almost  three  to  one. 

The  Nested  Integrals  method  is  more  efficient  than 
either  the  Rectangular  or  the  Trapezoidal  method.  The 
Nested  Integrals  results  were  as  accurate  as  the  Trapezoidal 
results  and  substantially  superior  to  the  Rectangular 
results.  Based  on  these  accuracy  and  efficiency  compari¬ 
sons,  the  Nested  Integrals  method  was  selected  for  further 
development  and  study. 
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A  trapezoidal  integration  approximation  rule  limits 
the  accuracy  of  both  the  Nested  Integrals  and  Trapezoidal 
algorithms.  An  accuracy  improvement  is  anticipated  in 
either  case  if  the  integration  is  approximated  using 
Simpson's  rule.  CRef  18:93-9*3  The  classical  debate  of 
whether  to  decrease  step  size  or  increase  integration  order 
can  be  raised.  For  the  problem  considered  here  and  prob¬ 
ably  for  other  more  complex  navigation  scenarios,  the 
sampling  theorem  compels  the  selection  of  step  size. 

Even  though  a  substantial  improvement  is  obtained  when  the 
rectangular  rule  is  replaced  by  the  trapezoidal,  the 
trapezoidal  algorithm  leaves  little  error  to  improve  on 
when  integration  step  size  is  well  below  the  Shannon  rate. 
For  this  reason,  an  extension  to  Simpson's  rule  for  the 
D(t)  computation  in  the  Nested  Integrals  algorithm  was 
not  made.  The  trapezoidal  rule  D(t)  approximation  is  the 
basis  for  the  Nested  Integrals  results  presented  in  the 
following  sections. 

3.6  Damped  Schuler  Case 

The  previous  undamped  Schuler  loop  example,  however, 
leaves  something  to  be  desired  in  terms  of  giving  insight 
into  practical  problems.  In  this  test  case,  the  variance 
being  approximated  grew  with  time;  therefore,  the  absolute 
error  could  also  grow  with  time  while  the  percent  error 
fell.  In  a  practical  case,  the  Schuler  loop  would 


probably  be  damped  by  measurements.  The  question,  then, 
is  whether  or  not  the  Nested  Integrals  method  still  gives 
valid  results  when  the  system  is  damped. 

To  answer  this  question,  two  chages  are  made  in  the 
problem  formulation.  First,  a  damping  term  corresponding 
to  continuous  velocity  measurement  and  feedback  is  added 
to  the  Schuler  Loop  model 


where  £  is  the  damping  factor.  For  this  numerical  example, 

C  is  set  to  0.3,  which  corresponds  to  an  underdamped  system. 
The  F  matrix 


F  = 


"  0  1  * 


yields  a  state  transition  matrix 


-CwJ-t-p) 

$(t,p)=«? 


{-(t)sv  sinwt) 


where 


(100) 


(cosoot+Cv  sinojt)  (-  —  sinwt) 


(cosujt-  £v  sincot) 

(101) 


to  =  to 


V  l-c2 


s 


V  =  1/V  1-  i* 


Second  the  gravity  disturbance  model  is  increased  in 
complexity  by  forming  a  second-order  filter.  This  gravity 
disturbance  yeilds  a  more  realistic  correlation  function. 
The  two-state  correlation  is  rounded,  as  empirical  data 
reflects,  at  zero  shift  rather  than  peaked  at  zero  shift 
as  is  the  correlation  for  the  previous  one-state  model. 

The  second  order,  stochastic  model  for  gravity  disturbance 


gives  an  autocorrelation  function 

(Xg-Xfl 

^>[6g(x1)6g(x2)]  =  a2«  d  (l  +  (102) 


The  parameter  d  is  no  longer  the  distance  at  which  corre¬ 
lation  is  l/e  of  the  zero  shift  value.  The  term  "corre¬ 
lation  distance"  of  Table  II  should  be  loosely  interpreted 
for  this  case.  Applying  (82)  and  (86)  yields 

?  -3(t9-t,|  /  \ 

QCr^hrU,,)]  =  o' e  1  l  1+(3|  t^t-J  J  (103) 

Constant  velocity  is  again  used  to  convert  the  spatial 
correlation  model  to  temporal; 


Let 


Pg(t) 


6g2  &g  y 

'  o2 

3a2* 

.yfig  y2  . 

.go2 

23 2a2_ 

(  104) 


for  filter  at  steady  state. 

An  augmented  system  can  again  be  formed  to  yield  a 

system  for  which  a  closed-form  solution  is  possible.  Let 
r 

6x 


Then 


= 


&g 
y  ^ 


(105) 


0 

1 

1  0 

1 

o' 

2 

~(l)c 

s 

-24^ 

1 

|  1 

0 

0 

0 

!  -3 

1 

0 

0 

!  0 

-3. 

*a(t)  +  «<— 


>  q(t) 


i 

v.  y 


(10  6) 


Let  F&  and  G&  be  defined  from  ( 106) ,  then  ( 14)  is  again 


valid.  Now, 


Qg  =  [43  V  ] 


(107) 


and  from  (67) 


Pa<0)=0‘ 


0 

0 

0 

0 


0 

0 

1 

3 


0 
0 
3 

23*J 


(108) 


For  numerical  results,  the  data  of  Table  II  was  again  used 
along  with  the  0.3  value  for  4  . 

Theoretically,  (14)  can  be  solved  in  closed  form 
using  (106),  (lO?),  and  (108).  For  this  example,  however, 
the  true  solution  was  found  by  solving  (14)  using  a  predictor- 
corrector  integration  algorithm  [Ref  14]. 

Using  (100),  (101),  (103),  and  (75) »  "the  Nested 
Integrals  results  for  position  and  velocity  variance  were 
produced  (Figures  D-l  through  D-4).  Figures  for  the  Damped 
Schuler  case  are  located  in  Appendix  D.  Since  these 
approximation  results  were  generally  more  accurate  than  for 
the  undamped  case,  an  additional  integration  step  size 
case  with  At  of  240  seconds  (6  on  the  logarithmic  scale) 
was  also  included.  The  approximation  algorithm  apparently 
does  solve  (29)  for  the  damped  Schuler  case  in  a  satisfac¬ 
tory  manner.  The  percent  position  variance  error  (Figure 
D-3)  would  be  acceptable  for  most  practical  situations, 
and  the  fact  that  this  error  stabilizes  and  does  not  grow 
is  reassuring.  The  percentage  results  for  the  undamped 
Schuler  case,  therefore,  were  not  dependent  on  the  pres¬ 
ence  of  a  ramp  term  in  the  variance  solutions. 

3. 7  Conclusions  from  Schuler  Loop  Cases 

Three  alternate  algorithms  were  produced  which  demon¬ 
strated  the  ability  to  solve  the  covariance  integral.  Each 
method  had  acceptable  small  error  in  calculating  P(t) 


84 


given  small  enough  integration  step  size  At.  The  Trapezoidal 
algorithm  was  more  accurate  than  the  Rectangular  algorithm. 
The  improvement  was  not  sufficient  to  warrant  developing 
even  higher  order  integration  algorithms.  On  the  other 
hand,  the  Nested  Integrals  algorithm  is  superior  to 
Trapezoidal  in  computation  efficiency  and  roughly  equiva¬ 
lent  in  accuracy.  The  Nested  Integrals  approach  was 
selected  for  further  development  and  study. 

The  undamped  Schuler  loop  used  in  this  method's 
comparison  left  some  doubt.  A  damped  Schuler  loop  case 
was  conducted  with  Nested  Integrals  to  demonstrate  that 
the  results  did  not  depend  on  the  undamped  model  insta¬ 
bility.  The  results  in  both  of  these  Schuler  loop  cases 
verify  that  the  Nested  Integrals  concept  is  a  valid  means 
of  producing  navigation  error  covariances. 

The  total  navaigation  system  error  propagation  in¬ 
cludes  more  than  p  Schuler  loop,  and  the  gravity  disturb¬ 
ance  is  also  multidimensional .  A  more  realistic  verifi¬ 
cation  is  required,  then,  to  demonstrate  the  performance 
with  a  complete  navigation  system  model. 


IV.  Nested  Integrals  Versus  Linear  State  Space 
Covariance  Analysis 

In  the  last  two  sections,  a  new  covariance  analysis 
algorithm,  Nested  Integrals,  was  developed  and  demonstrated. 
The  demonstrations  of  Section  III  were  one-dimensional. 

In  this  section,  an  analysis  of  a  complete  navigation 
system  is  presented  to  verify  further  the  validity  of  this 
method.  As  in  the  Schuler  loop  cases,  the  demonstration 
will  be  compared  to  a  linear  state  space  covariance  result. 
Since  Nested  Integrals  represents  an  alternative  to  the 
linear  state  space  covariance  analysis,  it  is  instructive 
to  demonstrate  not  only  where  the  methods  agree,  b”t  also 
where  they  yield  different  answers.  This  section  presents 
both  the  full-scale  verification  and  the  point  of  departure 
by  comparing  Nested  Integrals  results  to  those  produced 
by  linear  state  space  methods  on  a  range  of  trajectories. 

The  full-scale  verification  is  performed  using  a 
great  circle  trajectory.  On  this  trajectory,  the  linear 
state  space  covariance  analysis  theoretically  yields  the 
correct  result,  and  Nested  Integrals  answers  are  compared 
to  those  results  as  furtheb  proof  of  this  new  method.  Two 
minor  circle  trajectories  iare,  next,  selected  to  violate 
the  linear  state  space  method’s  restrictions.  The  Nested 
Integrals  method  properly  accounts  for  correlations  in  these 
minor  circle  instances,  and  these  results  should,  therefore, 
be  valid.  So  the  difference  between  the  two  methods  will 


be  the  error  induced  in  linear  state  space  covariance 
analysis  on  the  non-great  circle  trajectories.  This  para¬ 
metric  study  gives  insight  into  the  extent  of  mission 
trajectories  for  which  linear  state  space  results  might 
be  acceptable. 

To  conduct  these  studies,  models  for  the  trajectories 
mentioned  must  be  formed.  Also,  a  full-scale  gravity  dis¬ 
turbance  model  which  can  be  cast  in  a  linear  state  space 
format  must  be  defined.  With  a  definition  of  the  naviga¬ 
tion  error  propagation  model,  then,  the  study  can  begin. 

4. 1  Modeling  Choices 

The  model  requirements  for  both  types  of  analysis  are 
similar.  The  form  in  which  the  gravity-disturbance  statis¬ 
tical  model  enters  the  analysis  is  the  principal  difference. 
The  details  of  the  navigation  error  propagation,  gravity 
disturbance,  and  trajectory  models  are  provided  in  Appendices 
E  through  H.  The  following  brief  descriptions  provide 
additional  background  and  model  interface  detail. 

These  details  come  into  focus  when  the  needs  of  the 
two  analyses  are  reviewed.  Looking  first  at  the  linear 
state  space  covariance  analysis,  recall  the  structure  of 
this  method  from  Section  II i 

The  error  propagation  model  is  in  the  form 

x(t)  =  F(t)x(t)  +  G(t)u(t)  (3  ) 

The  gravity  disturbance  is  modeled  as  a  Gauss-Markov 


process  by 


V 


X  (t)  =  F  (t)x  (t)  +  G  (t)q(.t) 


(4  ) 


;g'  '*  g'  '-g  g 

The  disturbance  input  to  the  navigation  errors  is 

u(t)  =  Cg(t)Xg(t)  (  5  ) 

subject  to, 

=  Qg(t)6(t-p)  (6) 

The  augmented  state  is  defined  by 
x(t) 


= 


xg(t) 


so. 


xa(t) 


=  Fa(t)xa(t)  +  Ga(t)q(t) 


a'  _ 


where 


Fa(t) 


and 


'  F(t) 

G(t)  C^t) 

0 

F  ( t) 

g 

0 

(  7  ) 


(  8) 


(9) 


L  °g(t) 

The  covariance  of  the  augmented  state  satisfies 
Pa(t}  =  Fa(t)Pa(t)+Pa(t)F^(t)+Ga(t)Qg(t)GT(t) 
Solving  (14)  yields  the  linear  state  space  results 


(10) 


(14) 


The  state  x  and  x  must  be  defined  and  then  the  com- 

O 

ponents  of  each  associated  model  given.  The  navigation 
error  propagation  model  of  ( 3  ) »  for  this  study,  is  that 
of  Widnall  and  Grundy  [Ref  9  s 26-27],  The  vertical  channel 
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4'  -V. 


is  modified  to  represent  a  stable  system  based  on  altimeter 
aiding  as  suggested  by  Britting  [Ref  8  ].  The  navigation 
error  state  vector  of  the  Widnall-Grundy  model  has  nine 
elements 


x(t) 


(109) 


x9(t) 


which  are  defied  as  follows: 


X-^  is  6X,  error  in  longitude 
Xg  is  60,  error  in  latitude 
x^  is  6h,  error  in  altitude 

x^  is  6vg,  error  in  east  earth-relative  velocity 

x^  is  6vn,  error  in  north  velocity 

x^  is  6vz,  error  in  vertical  velocity 

Xr,  is  €  ,  orientation  error  about  the  east  axis 

i  G 

(local  level) 

x8  en’  or‘leni;a‘tlon  error  about  the  north  axis 

(local  level) 

x^  is  €z,  orientation  error  about  the  vertical  axis. 
The  modified  Widnall-Grundy  F-matrix  for  ( 3  )  is 
provided  in  Appendix  E.  Consistent  with  the  east-north- 
vertical  coordinate  frame  of  the  navigation  error  states, 
the  gravity  disturbance  is  defined 


f 

u(t)  =  6_gn(t)  =  6gn(t)  f 

K<*>  J 


(110) 
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where  the  n  superscript  indicates  the  physical  vector  is 
mathematically  represented  in  n-frame  ( east-north-vertical 
or  "navigation"  frame)  coordinates.  These  6gn  terms  drive 
the  respective  6vn  terms  in  ( 3  ) • 


G  = 


0 

i 

6 


(in) 


where  the  partitions  are  each  3^3  matrices. 

Since  a  linear  state  space  covariance  analysis  is  to 
be  performed,  a  statistical  model  for  gravity  disturbances 
in  the  form  of  ( 4  )  must  be  provided.  To  complement  the 
full-scale  navigation  error  model,  an  eight-state  gravity 
disturbance  model  is  selected  [[Ref  193 •  This  linear 
shaping  filter  has  been  the  subject  of  much  research  [Ref 
12l1  aimed  at  producing  a  model  which  replicates  empirically 
derived  anomaly  correlations  and  which  yields  auto-  and 
cross-correlations  of  other  disturbance  terms  consistent 
with  gravitational  field  theory.  The  filter  can  be  viewed 
as  three  separate  linear  systems  each  driven  by  an  element 
of  q.  Stationary  statistics  were  assumed  in  the  model 
development  [Ref  19 3,  so  the  noise  strength  Q  is  constant. 

D 

A  linear  combination  of  the  filter  states,  x  elements, 

O 

form  gravimetric  deflections  and  anomaly.  The  details  of 

the  C  (t)  output  matrix  of  (5  )  along  with  the  F  and  G 
g  g  6 

matrices  of  ( 4  )  are  provided  in  Appendix  F.  Q  is  also 

o 


described  in  Appendix  F. 


With  these  elements  of  the  analysis  defined,  the  F. 

cl 

and  G_  can  he  formed.  Thus,  P_  can  he  evaluated  hy  (14) 

CL  a 

and  numerical  solutions  for  P„  attained. 

a 

The  initial  condition  must  he  supplied  to  start  the 
solution  process.  Even  though  the  navigation  error  covari 
ance  is  zero  initially,  the  augmented  state  covariance  is 
non-zero,  because  the  modeled  gravity  process  is  assumed 
to  he  stationary,  hence  must  he  initialized  in  its  steady 
state  condition.  To  model  this  in-process  situation,  the 
gravity  disturbance  covariance  is  assumed  to  he  at  steady 
state  P  .  This  initial  condition  leads  to 


w  - 


o  :  o 


o  :  p. 


(112 


where  the  zero  matrices  are  of  required  dimension.  The  P 
initial  condition  for  the  gravity  disturbance  state's 
covariance  is  specified  in  Appendix  F. 

Now,  turning  to  the  Nested  Integrals  analysis,  the 
pertinent  equations  from  Section  II  ares 

x(t)  =  F(t)x(t)  +  G(t)u(t)  (3 


From  ( 28)  » 

Q(p.q)  =  ^Cu(p)uT(q) } 

From  ( 16) , 

^(t,^)  =  F(t)$(t,ti) 

for  t  >  t. ,  and  with  initial  condition  $(t, ,  t- )  =  I. 


(28a 


(16a 


(29a) 

(29b) 


P(t)  =  G(t)D(t)  +  DT(t)GT(t)+P(t)FT(t)+P(t)P(t) 

D(t)  =  S+  Q(t,p)GT(p)$T(t,p)dp 
Solving  (29)  yields  the  Nested  Integrals  result. 

F(t),  G(t),  and  u(t)  are  described  above.  The  Q(p,q) 
correlation  matrix  function  is  the  form  of  the  gravity 
disturbance  statistical  model  for  Nested  Integrals. 

Q(p,q)  is  formed  from  the  same  model  used  in  linear 
state  space  analysis.  Equation  (4  )  for  the  linear  state 
space  method  is  the  temporal  form  of  a  spatial  statistical 
process  (see  Appendix  B) .  The  correlations  of  Xg  are 
solved  in  this  spatial  domain  and  the  output  matrix  of 
(5)  applied  to  yield  the  elements  of  Q(p,q).  This  matrix 
function  is  provided  in  Appendix  F  with  the  method  for 
producing  the  central  angle  and  heading  associated  with  the 
r(p)  -  r(q)  pair. 

Solution  of  ( 29)  can  commence  with  the  specification 
of  zero  initial  condition  for  P(t)  as  discussed  in  Section 

II. 

The  linear  state  space  and  Nested  Integrals  analyses 
described  above  need  a  trajectory  model  to  evaluate  F(t), 
cg(t),  Qg,  Pg,tFg,  and  Q(p,q).  The  great  circle  trajectory 
is  described  in  Appendix  G  and  the  minor  circle  trajectory 
in  Appendix  H.  A  constant  altitude  profile  is  modeled  for 
both  the  great  circle  and  minor  circle  trajectories.  This 
flight  path  avoids  upward  continuing  [Ref  13]  the  gravity 
statistics.  For  convenience,  the  velocity  is  also  held 


constant  at  a  value  consistent  with  a  bomber  or  transport 
aircraft  cruise.  The  algorithm  generating  the  great  circle 
flight  path  (Appendix  G)  differs  from  the  minor  circle 
(Appendix  H)  in  that  an  inclination  can  be  specified.  The 
inclination  given  in  Table  IV  forces  the  great  circle 
path  off  a  cardinal  direction  and,  thus,  insures  that  all 
navigation  channels  are  exercised  in  this  study.  The  flights 
are  all  modeled  to  occur  in  the  same  geographic  area,  so 
minor  circle  results  can  be  compared  directly  to  great 
circle  results.  With  these  trajectory  models,  the  modeling 
structures  for  both  the  great  circle  and  the  minor  circle 
cases  are  set,  and  the  studies  themselves  can  be  addressed. 

4.2  Great  Circle  Case 

The  data  required  by  the  trajectory  model  of  Appendix 
G  and  the  gravity  disturbance  statistical  model  for  the 
great  circle  case  are  given  in  Table  IV  ,  below.  These 
data  will  be  used  in  several  studies  in  Sections  VI  through 
IX.  Also,  the  minor  circle  trajectories  will  use  the  same 
data  set  with  two  exceptions;  inclination,  i,  is  unused 
and  the  initial  circle  angle,  AQ,  has  a  different  meaning 
which  is  explained  in  Appendices  G  and  H.  The  Nested 
Integrals  integration  step  size  is  also  provided  in  Table 
TV. 

As  proposed,  both  Nested  Integrals  and  linear  state 
space  covariance  analyses  are  performed  on  the  great  circle 
case.  Since  the  trajectory  does  not  violate  the  linear 
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TABLE  LV  -  Great  Circle  Case  Data 


Quantity 

Symbol 

Value 

Units 

Final  time 

12240 

seconds 

Initial  time 

to 

0 

seconds 

Initial  circle  angle^ 

Ac 

0 

degrees 

1  _ 

Inclination  angle  * 

i 

^5 

degrees 

Velocity  magnitude* 

V 

615 

ft/sec 

Earth  radius 

Re 

20925640 

ft 

Altitude 

h 

0 

ft 

Gravity  magnitude 

g 

32.174 

ft/sec2 

Earth  rotation  rate 

^ie 

7.2921*10~5 

rad/sec 

Anomaly  variance 

<4 

1800 

(mgal)  2 

Correlation  parameter 

d 

20 

n  m 

Nested  Integrals  step  size 

At 

30 

sec 

f  Meaning  varies.  See  Appendix  G  for  great  circle  and 
Appendix  H  for  minor  circle  interpretations. 

*  Earth  relative  terms 

t  Applies  to  great  circle  only. 


state  space  constraints,  those  results  should  be  correct 
and  offer  a  benchmark  against  which  Nested  Integrals  per¬ 
formance  can  be  judged  on  this  full-scale  case.  The 
covariance  results  by  both  methods  at  200  minutes  into  the 
mission  are  provided  in  Table  V.  All  45  independent 
elements  of  the  covariance  matrix  are  given  in  lower  tri- 
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angular  format.  Using  the  linear  state  space  answers  as 
a  base,  the  differences  in  the  results  are  expressed  in 
percent  in  Table  V.  The  largest  difference  between  the 
two  methods  at  200  minutes  is  only  one-half  of  one 
percent. 

Table  V  gives  assurance  that  all  elements  are  being 
calculated  correctly  at  one  time.  Next,  circular-error- 
probable  is  calculated  throughout  the  mission  to  demonstrate 
that  the  solutions  agree  over  all  time.  The  circular- 
error-probable  (CEP)  summarizes  the  horizontal  position 
covariances  into  one  well-accepted  figure  of  merit.  This 
statistic  is  calculated  by  an  approximation  for  a  multi¬ 
dimensional  normal  distribution  [Ref  20  ];  xhe  calculation 
details  are  presented  in  Appendix  I. 

The  circular-error-probable  results  for  the  great 
circle  case  are  displayed  in  Figure  10.  Agreement  between 
the  two  analysis  methods  is  so  close  that  the  curve  markers 
must  be  offset  to  clarify  that  two  curves  are  plotted. 

These  results,  together  with  Table  V,  clearly  demonstrate 
the  validity  of  the  Nested  Integrals  analysis  on  a  full- 
scale  problem. 

In  this  case,  the  results  are  equivalent,  but  Nested 
Integrals  required  approximately  three  times  the  computa¬ 
tion  time  as  did  linear  state  space.  The  input/output 
time  for  Nested  Integrals  was  nearly  twice  the  computation 
time  because  of  the  21,012  file  GET  operations;  whereas, 
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Figure  10.  Circular-Error-Probable:  Great  Circle  Case 

Comparison 


the  linear  state  space  input/output  time  is  virtually  nil. 
While  no  efforts  were  made  to  minimize  computer  resource 
costs,  it  is  clear  that  a  substantial  penalty  is  paid  when 
Nested  Integrals  is  used  in  place  of  the  linear  state  space 
covariance  analysis;  therefore,  the  errors  produced  by 
linear  state  space  analysis  on  other  than  great  circle 
trajectories  must  be  considered. 

4.3  Minor  Circle  Cases 

Precisely  this  consideration  prompts  the  minor  circle 
trajectory  study.  The  great  circle  case,  viewed  as  a  3444 
n  m  minor  circle  radius  case,  has  shewn  exact  comparisons 
between  Nested  Integrals  and  linear  state  space  covariance 
analyses.  Two  other  minor  circle  cases  are  proposed  here 
to  scan  parametrically  the  trajectory  restriction. effects . 
Minor  circle  cases  are  defined  by  the  radius  of  the  closed 
flight  path  in  a  planar  sense.  The  details  of  the  minor 
circle  trajectory  are  provided  in  Appendix  H.  The  para¬ 
metric  range  should  include  trajectories  for  which  linear 
state  space  analysis  is  expected  to  err  significantly. 

Intuitively,  one  expects  little  error  on  cases  where 
the  minor  circle  radius  is  much  larger  than  the  correlation 
parameter.  Of  course,  any  trajectory  that  completes  a 
minor  circle  is  expected  to  elicit  some  error  when  using 
the  linear  state  space  method  since  the  linear  state  space 
modeling  grossly  misrepresents  the  correlations  along  such 
a  path.  More  pronounced  effects  should  be  observed  on 
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cases  where  the  minor  circle  radius  is  of  the  same  order 
of  magnitude  as  the  correlation  parameter.  To  cover  this 
range  of  effects,  additional  minor  circle  cases  of  19?  ran 
and  20  nm  radii  are  studied.  The  197  nm  case  completes 
exactly  one  circuit  in  the  204  minutes  of  simulated  flight. 
This  radius  is  approximately  five  times  the  20  ran  correla¬ 
tion  parameter,  so  this  case  should  represent  a  transition 
to  trajectories  with  significant  errors.  To  evoke  some 
definite  errors,  the  final  trajectory  is  a  relatively 
tight  turn  of  20  nm  radius.  This  radius  equals  the  corre¬ 
lation  parameter,  and  this  case  completes  nearly  ten  cir¬ 
cuits  during  the  204-minute  flight. 

The  parameteric  scan  of  this  study,  then,  includes 
trajectories  of  3444,  of  197,  and  of  20  nm  minor  circle 
radii.  Before  reviewing  the  covariance  results,  consider 
the  effect  this  set  of  trajectories  has  on  central  angle 
and,  therefore,  on  correlations.  The  central  angle  ijj ( t , tQ ) 
between  r(tQ)  and  r(t)  will  serve  as  the  example.  For  the 
minor  circle  cases  (see  Appendix  H  for  derivation), 


i|»(t,t  )  =  cos 


•1. 


Rg  -  rg[l-cos{(t-to)V/rc}] 


} 


(113) 


where  rc  is  minor  circle  radius  as  defined  in  Appendix  H. 

From  (113),  it  follows  that  the  great  circle  case,  r  =R  , 

C  6 

yields 

i|»(t,tQ)  =  |t-tQ|  V/rc  (114) 
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subject  to  0  <  \|»  <  180°.  This  central  angle  to  the  origi¬ 
nal  point  is  presented,  for  all  three  cases,  in  Figure  11. 
With  this  example  it  is  clear  that  only  the  great  circle 
case  gives  the  same  results  for  ^(tg.t^)  +  $('*'!» ^q)  as 
for  i|Kt2»t0)  —  a  matter  discussed  at  the  end  of  Section 
II  regarding  the  trajectory  restriction. 

This  spatial  variation  is  reflected  in  the  correla¬ 
tions  of  anomalous  gravity  terms.  For  illustration,  con¬ 
sider  the  correlation  of  anomaly  at  t  with  anomaly  at  t 
throughout  the  flights.  The  model  anomaly  correlation 
function  for  this  study  specifies  that  [Ref  19] 

tfggW  =  <Jg  [1  +  M  -  iM2]  e'M  (115) 

where 

M  »  r  ijj/d  (116) 

and 

r  is  position  vector  radius, 
d  is  correlation  parameter,  and 
o2  is  anomaly  variance. 

o 

All  three  quantities  are  defined  and  background  on  (115) 
is  given  in  Appendix  F.  The  correlation  of  Ag(tQ)  with 
Ag(t)  for  all  three  minor  circle  cases  are  presented  in 
Figure  12.  The  strong  effects  of  trajectory  on  correla¬ 
tion  are  clear.  On  every  complete  circuit,  the  anomaly 
correlation  with  the  initial  point  becomes  the  anomaly 
variance  again.  The  linear  state  space  covariance  analysi 
forces  the  great  circle  correlation  rule  regardless  of 


Figure  11.  Minor  Circle  Cases:  Central  Angle  Comparison 
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trajectory  and,  hence,  introduces  error. 

With  this  understanding  of  the  underlying  correla¬ 
tion  mismodeling,  the  results  should  come  as  no  surprise. 

The  3444  nm  case  results  were  presented  previously  in 
Figure  10.  The  197  nm  case  results  are  presented  in  Figure 
13,  and  the  20  nm  results  are  shown  in  Figure  14.  The  197 
nm  case  shows  a  close  agreement  between  the  two  analysis 
methods  over  the  first  half  of  the  flight,  as  one  would 
expect  from  the  apparent  match  between  the  great  circle 
and  the  197  nm  minor  circle  correlations  shown  in  Figure  12. 
The  correlation  curves  diverge  when  the  minor  circle  closes, 
and  a  slight  divergence  is  seen  in  the  circular-error- 
probable  results  toward  the  end  of  the  204-minute  mission. 
The  results  by  linear  state  space  on  the  197  nm  case 
appear  close  enough  to  the  current  Nested  Integrals  results 
for  practical  applications.  The  substantial  computational 
advantage  of  linear  state  space  would  motivate  the  accept¬ 
ance  of  such  minor  errors  on  this  mission. 

The  20  nm  circle  results  match  well  for  one  half  of 
the  minor  circle,  again.  In  this  case,  however,  the 
correlation  and  circular-error-probable  agreements  lasts 
only  one-twentieth  the  flight  time.  After  that,  the  linear 
state  space  covariance  analysis  result  is  drastically  dif¬ 
ferent  from  the  correct  Nested  Integrals  result.  The  two 
results  would  give  close  to  the  same  average  circular- 
error-probable  }  however,  the  peak  values  are  substantially 
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Figure  14.  Circular-Error-Probable t  20  nm  Minor  Circle 

Case 
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different.  While  the  linear  state  space  results  for  the 
197  nm  case  may  be  acceptable  as  an  approximation,  the 
results  on  the  20  nm  case  seem  hopelessly  wrong.  Each 
analyst  would  find  his  subjective  threshold  of  tolerable 
error  somewhere  in  the  span  of  minor  circle  radius.  The 
Nested  Integrals  approach  makes  visible  now  the  error  of 
the  linear  state  space  covariance  analysis. 

Another  insight  afforded  by  Nested  Integrals  analysis 
is  the  effect  of  the  trajectory  on  navigation  errors  pre¬ 
sented  in  Figure  15.  The  Nested  Integrals  results  are 
assumed  to  be  correct  due  to  the  great  circle  verification. 
A  Monte  Carlo  verification  of  these  minor  circle  results 
was  not  performed  since  a  more  useful  Monte  Carlo  verifi¬ 
cation  is  given  in  Section  V.  The  three  Nested  Integrals 
results  for  the  minor  circle  cases  are  plotted  on  the  same 
scale  for  comparison.  Note  that  the  higher  correlations 
on  the  smaller  minor  circles  do  not  breed  higher  errors  in 
general.  In  the  long  run,  the  great  circle  special  case 
produces  significantly  higher  errors. 


4.4  Comparison  Conclusions 

The  minor  circle  and  great  circle  studies  have  proved 
two  points  and  have  provided  some  insights  not  previously 
available.  First,  the  Nested  Integrals  algorithm  is  veri¬ 
fied  on  a  full-scale  problem.  Then,  the  potential  error 


of  using  linear  state  space  covariance  analysis  is  demon¬ 
strated  through  trajectory  variations  out  of  the  great 
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Figure  15.  Minor  Circle  Cases  Comparison,  Nested 

Integral  Results 


circle  restriction.  This  latter  process  gives  insight  into 
how  significant  these  errors  might  he  and  into  when  the 
error  level  might  he  accepted  in  order  to  gain  the  greater 
computational  efficiency  of  linear  state  space  methods. 
Finally,  a  comparison  of  gravity-induced  navigation  errors 
as  a  function  of  the  trajectory  is  provided  by  the  Nested 
Integrals  method.  The  results  are  interesting  and  provide 
a  view  of  issues  not  seen  before.  The  fact  that  Nested 
Integrals  can  correctly  calculate  these  covariances  in  an 
efficient  manner  is  a  prime  aim  of  this  research.  A  needed 
alternative  to  the  linear  state  space  analysis  trajectory 
restrictions  has  been  presented  on  a  full-scale  problem. 


V.  Nested  Integrals  Versus  Monte  Carlo  Analysis 

The  three  statistical  analyses  discussed  in  Sections  I 
and  II  are  linear  state  space  covariance  analysis,  Monte 
Carlo  analysis,  and  covariance  integral  analysis.  Nested 
Integrals,  as  a  covariance  integral  type  of  analysis,  was 
compared  to  linear  state  space  analysis  in  Section  IV. 

The  minor  circle  trajectories  used  in  Section  IV  while 
violating  the  linear  state  space  restrictions  still  do  not 
exercise  the  full  range  of  Nested  Integrals  capability. 

Also,  the  question  of  whether  this  new  method  is  more  effi¬ 
cient  than  Monte  Carlo  must  be  faced.  For  these  reasons, 
this  section  is  devoted  to  a  Nested  Integrals  to  Monte 
Carlo  comparison  on  a  complex  trajectory.  The  vehicle  for 
this  comparison  is  a  recently  completed  Monte  Carlo  study 
of  an  air-launched  strategic  missile  mission  [Ref  10]. 

Nested  Integrals  method  is  used  on  the  same  problem  allowing 
a  direct  comparison  of  results  and  computational  costs. 

The  comparison  of  Nested  Integrals  to  Monte  Carlo 
method  is  the  primary  motivation  for  this  particular  study, 
but  several  other  features  are  offered: 

a.  Complex  trajectory.  The  strategic  missile  trajec¬ 
tory  is  significantly  different  from  the  cons tan t-groundspeed, 
level-flight  mission  used  in  Section  IV.  While  the  missile 
trajectory  stays  in  a  great  circle  plane  the  flight  profile 
has  significant  acceleration  and  altitude  changes. 

Groundspeed  increases  some  3285  ft/sec  while  altitude 
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changes  almost  80,000  feet  in  the  initial  boost  and  climb 
phases.  Again,  the  trajectory  flexibility  required  for 
studies  such  as  this  is  a  prime  motive  for  Nested  Integrals 
development. 

b.  Correlation  model  variation.  The  Monte  Carlo 
analysis  employed  a  different  statistical  model  [Ref  21] 
than  the  previous  linear  state  space  model.  For  this  com¬ 
plex  altitude  profile,  a  statistical  model  must,  and  this 
one  does,  give  correlations  at  and  between  different  alti¬ 
tude  levels  which  are  consistent  with  the  interrelation¬ 
ships  imposed  by  gravitational  field  theory  (upward  continu¬ 
ation)  .  The  flexibility  to  use  a  variety  of  statistical 
models  is  also  one  prime  motive  for  the  Nested  Integrals 
development. 

c.  Real-world  problem.  The  problem  addressed  here 
is  not  academic.  Analyses  of  this  type  are  typical  of 
the  exercise  performed  in  generating  an  error  budget  for  a 
proposed  new  system.  In  this  sense,  the  air-launched 
strategic  missile  represents  a  real-world  problem  on  which 
the  Nested  Integrals  methods  can  be  applied. 

d.  Independent  verification.  The  author  was  the 
Air  Force  study  manager  for  the  Monte  Carlo  analysis,  but 
the  model  selection  and  analysis  were  performed  by  the 
contractor,  a  separate  party.  These  published  results 
[Ref  lo]  offer  an  independent  check  on  the  newly  developed 
Nested  Integrals  method. 
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The  rationale  for  this  comparison  is  to  emulate  the 
Monte  Carlo  study  as  nearly  as  possible.  The  Monte  Carlo 
study  is,  therefore,  discussed  before  the  Nested  Integrals 
model  selection.  Finally,  navigation  accuracy  results 
are  compared  and  computational  costs  are  contrasted. 

5.1  Advanced  Strategic  Air-launched  Missile  Problem 

The  problem  is  the  medium  through  which  Nested  Integrals 
is  compared  to  Monte  Carlo  analysis.  The  genesis  of  the 
problem  was  a  system  accuracy,  error- budget  exercise  con¬ 
ducted  by  Aeronautical  Systems  Division  of  the  Air  Force 
Systems  Command.  Several  trajectories  were  considered 
in  an  elaborate  trade-off  of  systems  and  flight  strategies. 

The  linear  state  space  covariance  analysis  approach 
is  typically  used  for  estimating  the  gravity  model  contri¬ 
butions.  This  trajectory  clearly  violates  the  constant 
altitude  model  restriction,  and  some  analysis  considering 
upward  continuation  was  desired.  Geodynamics  Corporation 
was  tasked  in  August,  1977  to  perform  an  analysis  which 
would  give  the  gravity  model  contribution  to  inertial  navi¬ 
gation  errors  on  this  complex  mission.  The  results  of  that 
study  are  recorded  in  Reference  10  dated  August,  1978.  Of 
the  studies  included  in  that  work,  the  1500-nm  missile 
trajectory  is  selected  as  a  case  on  which  to  apply  Nested 
Integrals  analysis. 
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5.2  Monte  Carlo  Stud: 


The  aspects  of  the  Monte  Carlo  study  which  affect  the 
course  of  the  Nested  Integrals  analysis  are: 

a.  the  trajectory, 

b.  the  navigation  simulation,  and 

c.  the  statistical  model. 

These  categories,  obviously,  are  the  three  models  needed 
for  Nested  Integrals  analysis. 

The  trajectory  used  for  the  Monte  Carlo  study  is 
described  on  page  6  of  Reference  10.  The  necessary  excerpts 
from  that  description  are  presented  in  Table VI,  below. 

The  Monte  Carlo  trajectory  generator  used  acceleration 
polynomials  to  match  segment  endpoint  conditions  [[Ref  10:4]. 
The  position  time  history,  0Q,  was  calculated  from  the 
resulting  acceleration  profile.  With  both  acceleration  and 
position  profiles  for  the  design  mission,  the  specific 
force  profile  was  generated  using  a  truth  model  for 
gravity  [Ref  10:4] 

f(t)  =  a(t)  -  G[r(t)]  (  117) 

where  f(t)  is  specific  force, 

a(t)  is  total  acceleration, 

G[]  is  true  gravitation,  a  vector  not  to  be  confused 
with  the  distribution  matrix  of  the  linearized 
navigation  error  propagation  model,  and 

r(t)  is  from  6qI  the  design  mission. 

The  navigation  simulation  was  whole  valued  and  based  on 

solving 


Table  71  -  Missile  Flight  Path  Data 
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These  quantities  are  specific  inputs  to  PROFGEN  CRef  7]  and  are  discussed 
in  subsection  5-3- 


Letting  "*■"  represent  navigation  estimates,  the  navigation 
simulation  solved 

r(t)  =  f(t)  +  Gm[  r(t)]  (118) 

where  Gm  is  the  gravitational  model  simulated  here  as 

£m[r(t)]  =  G[r(  t)  ]  +  6g[r(t)]  (119) 

and  6g[r( t) ]  is  the  gravity  disturbance  stochastic  reali¬ 
zation  from  the  statistical  model.  The  navigation  position 
error  is  defined  as 

6r(t)  =  r(t)  -  r(t)  (120) 

The  Monte  Carlo  approach  is  to  produce  an  ensemble  of 
6r(t)'s  based  on  an  ensemble  of  6g[_r(t)]'s.  The  gravity 
disturbance  ensemble  is  generated  in  such  a  manner  that, 
in  the  limit,  the  correlations  of  the  statistical  model  are 
replicated.  The  theory  is  that  the  ensemble  of  naviga¬ 
tion  error  time  histories  is  representative  of  those  which 
would  be  produced  by  the  population  of  6g’s  whence  the 
original  statistics  were  derived.  For  this  study  Geodynamics 
defined 


6gf(ti) 


-< 


1 


(vAgi  +  2(g/r)NiJ 


(121)* 


*  This  notation  is  maintained  here  for  traceability  to  the 
Monte  Carlo  study.  Crossrange  deflection  m  is  equivalent 
to  transverse  deflection  p,  and  /  to  r  from  Appendix  F. 
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where 

/  is  the  downrange  deflection 

m  is  the  crossrange  deflection 

Ag  is  anomaly 

N  is  geoidal  height 

t  superscript  indicates  the  local-level  downrange- 
crossrange -vertical  coordinate  frame. 

The  2(g/r)N^  term  corrects  for  the  fact  that  Ag^  is 

not  exactly  6g„ .  This  minor  correction  has  an  insignif i- 

cant  effect  on  the  results  calculated  in  this  case  which 

is  why  6g  =  -Ag  was  used  as  the  model  in  Section  IV.  Due 
z 

to  computer  core  limitations,  the  desire  for  more  dense 
sampling  of  the  deflection  disturbances,  and  an  assumed 
altimeter-aided  system,  this  hypothetical  vertical  distur¬ 
bance  was  not  even  simulated  [Ref  10:493  in  the  Monte  Carlo 
study. 

Next,  a  statistical  model  was  needed  on  which  to 
base  the  ensemble  of  gravity  disturbance  profiles.  The 
Tscheming-Rapp  anomaly  degree  variance  model  £Ref  21:43- 
463  number  4  was  selected.  The  mathematical  details  of 
this  model  are  repeated,  in  part,  in  Appendix  J.  This 
statistical  model  is  based  on  a  spherical  harmonic  repre¬ 
sentation  of  the  anomaly  autocorrelation,  assumed  to  be 
spatially  stationary  and  isotropic. 

This  model's  coefficients  explicitly  represent  the 
contribution  to  anomaly  variance  of  spatial  frequency 
terms  -  hence  the  title  anomaly  degree  variance  model.  A 
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closed  mathematical  form  relating  harmonic  coefficient  to 
degree  number  is  hypothesized  (see  Appendix  J)  and  the 
resultant  model  fit  to  empirical  data  by  a  parameter 
identification  process.  Using  gravitational  field  theory 
the  auto-  and  cross-correlaticns  of  and  between  deflection 
and  geoidal  undulations  are  derived  in  turn  using  the 
anomaly  correlation  model  as  a  basis.  Throughout  this 
derivation,  the  altitude  coordinate  information  for  both 
positions  of  the  correlation  evaluation  is  maintained  and, 
thereby,  provides  the  theoretically  consistent  means  of 
upward  continuing  the  statistics. 

To  employ  this  model,  &g  evaluations  were  generated 
for  98  position  points  corresponding  to  25-second  time 
increments.  Let  i  subscripts  represent  t^  quantities  for 
i=l,  ....  98.  Then  define 

6Ii  =  6f>(Ei.)  (122) 

From  the  sequence  {Sg^}  =  [Sg-^,  6g2*  *•*  de;fine  "the 

ensemble  mathematical  vector 
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Let  0  represent  the  correlation 

T  =  &gT] 


(123) 


(124) 


Equation  (124)  can  be  evaluated  given  the  statistical 
model,  (122),  and  the  sequence  £r^}.  Define  a  square  root 
matrix  S  by 

T  =  SST  (  12  5) 

Then  a  population  having  the  same  covariance  0  can  be 
generated  by 

6g"  =  Sq  (126) 

where  comes  from  a  multidimensional  gaussian  population 
with  properties: 

$ [3 ,]  -  0  (127) 

and 

(gWl  -  1  C12S) 

An  ensemble  of  90  q  vectors  was  used  to  generate  an 
ensemble  of  gravity  disturbance  sequences  using  (126)  and 
(123).  The  disturbance1  is  introduced  through  (119)  to  the 
navigation  solution  of  (118).  The  result  is  an  ensemble 
of  navigation  error  time  histories,  each  given  by  (120), 

From  the  ensemble  of  6r(t),  the  horizontal  components 
were  processed  to  yield  downrange  and  crossrange  variances 
at  each  point.  Let  and  oQ  represent  these  variances. 
Then  circular-error- probable  (CEP)  was  calculated  from 


CEP  =0.562  omax  +0.615  Omin 

(129) 

where 

amax  ■  Maltimum  t°d-  ae] 

(130) 

°min  ■  Minimum  C°d"  oc} 

(131) 
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5*3  Trajectory  Model  for  Nested  Integrals 

With  the  Monte  Carlo  results  established,  attention 
can  be  turned  to  the  modeling  choices  for  Nested  Integrals. 
Obviously,  the  statistical  model  must  be  the  same  as  for 
Monte  Carlo.  The  trajectory  model  should  be  the  same,  but 
the  polynomial-fit  trajectory  generator  was  not  available 
for  the  Nested  Integrals  analysis.  The  trajectory  was 
simulated  using  PROFGEN  [Ref  ?  j,  a  computer  program 
designed  to  use  trajectory  segment  definitions  as  input 
and  give  a  complete  mission  time  history  as  output.  The 
PROFGEN  segment  descriptions  and  results  are  also  provided 
in  Table  VI  for  direct  comparison  to  the  Monte  Carlo  trajec¬ 
tory  data.  The  six  segments  are  as  follows: 

a.  Launch  and  accelerate  (5*137  g's)  from  615  to 
2600  ft/sec  groundspeed 

b.  ’  Pitch-up  vertical  turn  (7  g's)  to  a  pitch  angle 
of  16.3° 

c.  Accelerate  (5*137  g’s)  while  climbing  until  3900 
ft/sec  groundspeed. 

d.  Continue  straight  climb  at  constant  groundspeed. 

e.  Pitch-over  vertical  turn  (7  g's)  level  out  at 
80000  ft  altitude. 

f.  Cruise  at  constant  altitude  and  groundspeed. 

The  terminal  dive  phase  of  the  last  few  seconds  of  simu¬ 
lated  flight  was  not  modeled  for  the  Nested  Integrals  analy¬ 
sis. 
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5«4  Inertial  Navigation  Error  Propagation  Model  for 

Nested  Integrals _ 

The  inertial  navigation  simulation  for  Monte  Carlo 
was,  in  essence,  that  of  an  inertial  computation  frame. 

The  choice  of  instrument  or  computation  frame  does  not 
change  the  gravity  disturbance  vector  and  Britting  [Ref  8] 
has  shown  a  general  equivalence  of  error  propagation  models; 
therefore,  the  modified  Widnall-Grundy  model  used  in 
Section  IV  was  used  for  this  analysis  as  well.  The  verti¬ 
cal  channel  disturbance  was  eliminated,  and  the  Nested 
Int  cals  step  size  was  set  at  25  seconds  to  be  consistent 
with  the  modeling  and  sampling  of  the  Monte  Carlo  study. 
Circular- error-probable  was  calculated  from  the  Nested 
Integrals  covariance  results  using  the  same  downrange- 
crossrange  rule  as  for  Monte  Carlo. 

5 • 5  Comparison  of  Results 

The  circular-error-probable  results  for  both  the 
Monte  Carlo  and  the  Nested  Integrals  analyses  are  presented 
in  Figure  16.  The  95%  confidence  bounds  for  the  Monte  Carlo 
results  are  also  repeated  [Ref  10]  here. 

The  navigation  accuracies  calculated  by  these  two 
statistical  methods  differ  somewhat  due  to  modeling  differ¬ 
ences  and  Monte  Carlo  sample  size  limitations.  The  results 
are  generally  equivalent,  however,  and  another  degree  of 
confidence  in  the  Nested  Integrals  approach  is  warranted. 
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Figure  16.  Circular-Error-Probable t  Air  Launched 
Strategic  Missile  Case 
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Time  (min) 


The  computational  costs  were  quite  different.  For 
the  90-case  ensemble,  the  Monte  Carlo  computational  time 
was  30  times  the  Nested  Integrals  time.  That  is,,  the 
entire  Nested  Integrals  analysis  was  performed  in  the  same 
computation  time  required  for  three  of  the  Monte  Carlo 
samples.  Data  was  not  available  to  compare  post-processing 
time,  but  the  Monte  Carlo  study  requires  significantly 
more  post-processing  since  the  sample  outputs  are  6r(t),  not 
covariance  per  se. 

These  results  demonstrate  that  Nested  Integrals  has 
fulfilled  the  expectations  expressed  in  Section  II.  The 
statistical  analysis  is  as  flexible  as  the  Monte  Carlo 
method,  but  significantly  less  expensive  computationally. 

The  numerical  results  by  the  two  methods  are  equivalent;  so, 
the  Nested  Integrals'  computational  efficiency  does  not 
entail  any  loss  of  authority  in  the  results. 


123 


VI .  Input  Correlation  Function  Variations 

In  Sections  IV  and  V,  full-scale  navigation  error 
analyses  were  performed  using  Nested  Integrals  but  based 
on  different  gravity  disturbance  statistical  models. 

Nested  Integrals  method  is,  in  part,  motivated  by  the  need 
to  consider  different  statistical  model  forms.  This  sec¬ 
tion  gives  a  direct  comparison  of  results  from  these  two 
models  and  from  another  fully  developed  model  which  has 
not  been  previously  presented  in  this  work.  The  purpose 
of  this  comparison  is  two-fold: 

a.  First,  to  demonstrate  the  capability  to  perform 
analyses  with  a  variety  of  statistical  model 
forms . 

b.  Secondly,  the  linear  state  space  model  has  so 
dominated  past  studies  that  the  existence  of  other 
models  must  be  emphasized. 

This  simple  comparison  should  demonstrate  the  availability 
of  fully  developed  alternate  statistical  models. 

The  demonstration  requires  for  each  analysis  the 
usual  three  models:  trajectory,  error  propagation,  and 
disturbance  correlation.  The  disturbance  correlation  model, 
as  discussed,  will  be  varied  with  three  separate  mathe¬ 
matical  forms  considered.  These  three  statistical  models 
are  discussed  below  in  subsections  6.1,  6.2,  and  6.3.  For 
each  correlation  model,  a  complete  Nested  Integrals  analysis 
is  performed  based  on  the  same  trajectory  and  error  propa- 
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gation  models.  That  is,  everything  in  the  analysis  is 
held  to  a  constant  form  except  the  correlation  model.  The 
benign  great  circle  trajectory  of  Section  IV  is  used  along 
with  the  modified  Widnall-Grundy  error  propagation  model. 

The  necessary  data  for  these  two  models  are  found  in 
Table  IV  of  Section  IV* 

6.1  Linear  State  Space  Correlation  Model 

This  model  is  presented  in  Appendix  F,  and  its  use 
has  been  discussed  in  Section  IV.  Some  additional  comment 
is  due  since  this  model  is  being  compared  here  with  other 
model  forms. 

The  basis  of  this  modeling  technique  is  a  Gauss-Markov 
process.  The  gravity  disturbances  are  represented  as  the 
outputs  of  a  linear  filter  driven  by  white  gaussian  noises. 
The  resulting  disturbances  quantity  auto-  and  cross-corre¬ 
lations  are  consistent  with  gravitational  theory  except 
for  minor  approximations  of  the  correlations  between  anomaly 
and  undulation  and  between  anomaly  and  downrange  deflec¬ 
tion  [[Ref  195B-6].  The  model  used  in  Section  IV  and  here 
is  based  on  an  eight-state  filter  driven  by  three  indepen¬ 
dent  noises.  The  resulting  anomaly  autocorrelation  is 

«*gg(£.£')  =  0gg(ijO  =  Og  (1  +  M  -  iMZ)  e~M  (115) 

where 

M  =  r  v(j/ d  (116) 

and 
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4r  is  the  central  angle  separating  r  and  r*, 
r  is  the  reference  radius  value,  and 
d  is  the  correlation  parameter. 

As  it  stands,  this  model  is  only  valid  for  one  altitude, 
so  the  correlation  is  shown  as  a  function  of  central  angle 
shift  alone. 

o 

The  anomaly  variance  level  of  1800  mgal  used  for 
this  model  is  consistent  with  the  worldwide  anomaly  vari¬ 
ance  cited  in  References  22  and  23.  The  correlation  para¬ 
meter  of  20  nm  (Table  IV)  is  an  arbitrary  figure  selected 
for  the  Section  IV  study.  No  one  value  of  d  will  yield  an 
overall  good  fit  to  empirical  data.  Several  independent 

models  of  this  type  could  be  postulated,  each  with  d^  and 
2 

a  parameters.  The  linear  state  space  disturbance  model 
6i 

used  here  should  be  viewed  as  a  representative  example,  not 
an  ultimate  model.  The  empirical  correlation  function  of 
Reference  23  has  a  correlation  distance  of  41  nm,  so  the 
20-nm  correlation  parameter  of  this  model,  will  yield  dif¬ 
ferent  results  on  that  account. 

The  correlation  distance  of  the  anomaly  degree  vari¬ 
ance  model  is  closer  to  t]  * s  20  nm  parameter  but  the  vari¬ 
ance  level  is  lower.  The  navigation  error  correlation 
results  are  linear  with  variance  level,  assuming  the  entire 
anomaly  correlation  is  simultaneously  scaled  up  or  down 
with  the  variance  changes.  The  effect  of  correlation  para¬ 
meter  change  is  trajectory  dependent.  Generally  one  expects 
disturbances  from  a  longer  correlation  parameter  model  to 


appear  more  like  a  constant  input  to  the  navigation  algo¬ 
rithm  which  acts  as  a  low  pass  filter.  The  shorter  corre¬ 
lation  parameter  will  induce  higher  frequency  inputs  and 
tend  to  yield  attenuated  errors. 

The  linear  state  space  correlation  model  example  has 
the  same  trajectory,  error  propagation  and  disturbance 
statistical  models  as  the  great  circle  case  of  Section  IV. 
The  circular-error-probable  results  are  identical  to  the 
Section  IV  study  and  are  repeated  here  for  comparative 
purposes. 


6.2  Anomaly  Degree  Variance  Correlation  Model 

This  model  was  introduced  in  the  Monte  Carlo  study 
discussed  in  Section  V.  Further  details  on  this  modeling 
technique  are  presented  in  Appendix  J. 

The  anomaly  correlation  function  is  assumed  to  be 
stationary  and  isotropic.  The  function  can,  therefore, 
be  expanded  in  spherical  harmonics  (i.e.  Legendre  functions). 
The  symmetry  of  the  isotropic  correlation  functions  means 
the  harmonic  expansion  can  be  limited  to  the  zeroth  order 
(i.e.  Legendre  polynomials). 


w  /  r2  \n+2 

«gg(E'£'>  =  =  n?0  cn(^)  Pn(00S  ,f) 

(  132) 

The  anomaly  variance  is  given  by 

0*  =  2  c  (133) 

S  n=o  n 


128 


based  on  setting,  r=r'=R  and  ip=0 .  A  mathematical  expres¬ 
sion  is  hypothesized  for  the  manner  in  which  cn  approaches 
zero  as  n  approaches  infinity.  Each  cn  represents  the 
contribution  to  the  variance  from  the  n-th  degree  poly¬ 
nomial,  hence  the  name  anomaly  degree  variance  model.  Note 
that  the  spherical  harmonic  mathematical  form  inherently 
gives  upward  continuation  of  the  statistics  without  an 
additional  integration.  Closed-form  mathematical  expres¬ 
sions  for  ( I32)  are  developed  and  associated  disturbance 
quantity  correlations  derived  from  this  basis  (see  Appen¬ 
dix  A) . 

The  specific  cn  model  and  the  associated  parameter 

set  are  discussed  in  Appendix  J.  The  resulting  model  has 

2 

an  anomaly  variance  of  1795  mgal  and  a  correlation  dis¬ 
tance  of  approximately  30  nm.  The  circular- error-probable 
for  the  great  circle  trajectory  for  this  statistical  model 
is  also  presented  later. 

6.3  Attenuated  White  Noise  Correlation  Model 

The  two  previously  discussed  models  were  introduced 
prior  to  this  Section.  Another^fully  developed  statistical 
model  exists  and  is  included  in  this  model  comparison. 

The  attenuated  white  noise  model  is  presented  in 
Reference  23  with  additional  discussion  in  Reference  22. 
Some  of  the  details  are  reiterated  here  in  Appendix  K. 

The  basis  of  this  model  is  a  hypothetical  white  noise, 
anomalous  potential  process  on  a  spherical  shell  at  a  depth 
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d  below  the  earth’s  surface.  The  earth-surface  anomalous 
potential  is,  therefore,  the  upward  continuation  of  this 
white  noise  process.  LaPlace’ s  equation  relates  anomalous 
potential  above  the  model  noise  process  and  the  resultant 
earth-surface  process  is  correlated  by  way  of  the  Poisson 
Integral  upward  continuation  [Ref  3]* 

The  mathematical  expressions  which  result  from  this 
hypothesis  are  quite  complex  and  an  asymptotic  form  is 
postulated.  The  asymptotic  form  maintains  accuracy  within 
an  earth's  radius  of  the  surface,  and  this  region  includes 
all  "terrestrial"  navigation  missions.  The  anomaly  corre¬ 
lation  from  the  asymptotic  form  of  the  attenuated  white 
noise  model  is 

0gg(r.r’)  =  0gg(^»  h»  h') 

8d4(2d+h+h')[2(2d+h-t-h’)2  -  3(Rfr)2]  2  ,  . 

[(2d+h+h')2  +  (Rijj)2]7/2  °&  13 

where 

is  central  angle  between  r  and  r' , 
h,h'  are  geocentric  altitudes  from  r,  r'  respectively, 
d  is  white  noise  depth, 

2 

a  is  anomaly  variance,  and 

o 

R  is  radius  of  sphere  which  approximates  earth's 
surface . 

Note  that  the  model  has  upward  continuation  of  statistics 
by  the  altitude  data  in  the  functional  expression.  To  fit 
empirical  data,  three  statistically  independent  white  noise 
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shells  were  postulated  each  with  a  variance  level  and  depth. 
This  device  allows  both  long  and  short  wavelength  informa¬ 
tion  in  the  empirical  function  to  be  modeled.  The  final 

correlations  are  the  sum  of  these  three  independent  models. 

2 

With  the  parameter  sets  (d. }  and  {<j_  }  identified 

1  gi 

(Appendix  K) ,  the  resulting  model  has  an  anomaly  variance 
2 

of  1821  mgals  .  The  results  using  this  statistical  model 
on  the  great  circle  trajectory  are  also  presented  in 
Figure  17 . 

6.4  Comparison  of  Navigation  Errors 

The  three  model  concepts  presented  above  attempt  to 
replicate  empirical  statistics  and  to  conform  to  the  stric¬ 
tures  of  gravitational  field  theory.  The  basic  premise 
for  the  mathematical  form  is  different  in  each  case,  and 
the  result  is  different  representations  for  each  correla¬ 
tion  of  interest.  Equations  (115),  (132),  and  (134)  pre¬ 
sent  the  different  functional  forms  for  anomaly  autocorre¬ 
lation  as  examples.  Each  model  employs  the  isotropic  assump¬ 
tion,  so  only  four  disturbance  correlations  are  required. 
These  four  correlations  are  presented  in  Appendix  L  in  a 
comparative  study  of  the  three  models  above. 

The  circular-error-probable  results  for  each  model  on 
the  great  circle  trajectory  are  presented  in  Figure  17. 

These  results  are  considerably  different  as  the  Appendix  L 
comparison  forbodes.  Which  model  is  correct?  There  is  no 


general  answer  to  that  question.  Each  model  purports  to 
represent  the  same  physical  process,  and  each,  being 
limited,  fails  in  some  way.  Which  model  is  "correct"  is  a 
subjective  judgment  that  will  depend  on  the  mission  being 
analyzed  and  the  analyst's  prejudices. 

The  attenuated  white  noise  and  the  anomaly  degree 
variance  models  have  upward  continuation  of  the  statistics 
inherent  in  the  mathematical  form,  a  desirable  general 
feature.  These  two  models  are  also  internally  consistent 
with  gravitational  field  theory.  The  linear  state  space 
model  requires  additional  work  to  upward  continue  statis¬ 
tics  [Ref  13]  and  only  approximates  consistency  with  gravi¬ 
tational  field  theory. 

The  outcome  of  this  investigation  will  not  specify  a 
model  type.  This  study's  aims  are  simply  to 

a.  Point  out  the  existence  of  statistical  model 
alternatives,  and 

b.  To  demonstrate  Nested  Integrals  method's  capa¬ 
bility  to  employ  any  of  these  different  mathe¬ 
matical  forms. 

The  analysis  results  in  Figure  1 7  provide  the 
demonstration  of  Nested  Integrals  capability  and  simultane¬ 
ously  shows  that  model  selection  has  a  direct  effect  on 
analysis  results. 


VII .  Variations  in  Spherical  Harmonic  Modeling 

Inertial  navigation  gravity  model  improvements  should 
be  evaluated  in  terms  of  resulting  system  true  accuracy. 
Model  improvement  changes  the  residual  field  statistics 
both  in  magnitude  and  in  spectral  content.  If  these 
changes  can  be  summarized  in  residual  field  correlation 
models,  Nested  Integrals  offers  an  evaluation  technique  to 
select  the  model  complexity  which  meets  mission  accuracy 
specifications.  This  section  demonstrates  such  an  evalu¬ 
ation  on  a  hypothetically  spherical  harmonic  model  of 
limited  degree  and  order  and  with  no  error  in  harmonic 
coefficients.  While  such  a  model  is  unrealistic,  a  study 
based  on  this  hypothetical  model  can  give  insight  into  the 
sensitivity  of  system  accuracy  to  model  complexity. 

As  in  Section  VI,  this  analysis  will  involve  varia¬ 
tions  in  the  input  correlation  function.  Unlike  Section 
VI,  the  functional  form  will  be  constant.  For  this  study, 
the  model  changes  are  seen  in  model  parameter  variations. 

The  correlation  modeling  is  discussed  below.  The  great 
circle  trajectory  model  and  modified  Widnall -Grundy  error 
propagation  model  are  again  used  a  vehicle  for  the 
statistical  model  variations. 

7.1  Correlation  Model 

The  gravitational  model  can  be  based  on  a  spherical 
harmonic  expansion  of  gravitational  potential  U  CRef  3  s342] 


FAQ I  BUMUNOT  FUND 


135 


U(r) 


JnPn^C0S  ^ 
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j  =  integer 
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T(t)  =  P  f  t) 
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where 


|i@  is  earth  gravitational  constant 
r  is  radius  vector  of 

a  earth  semi-major  axis 

0  colatitude 
X  longitude 

Jn  are  zonal  harmonic  coefficients 

J _ ,K _  are  tesseral  harmonic  coefficients 

nra  nra 

Pn( • )  are  Legendre  polynomials  of  degree  n 

Pnm( • )  are  Legendre  functions  of  degree  n  and  order  m 

Zonal  harmonic  models  of  low  degree  and  order  are  most 
common,  but  higher  order  models  have  been  formed  CRef  103- 
Specific  values  for  the  parameters  depend  on  the  geodetic 
model  selected.  This  work  does  not  depend  on  any  specific 
model  so  these  data  are  not  required  explicitly  here.  Models 
are  based  on  truncated  versions  of  (135)  making  use  of 

G(r)  =  v  U(r)  (139) 
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where  G  is  the  gravitational  vector.  Any  attempt  to 
identify  coefficients  will  result  in  errors  from  the  fol¬ 
lowing  sources: 

a.  Measurements  are  limited  by  economics  and 
politics,  hence  aliasing  of  higher  order  effects 
onto  low  order  coefficients. 

b.  Measurement  errors. 

Koch  [Ref  24]  has  pointed  out  rather  severe  aliasing  at 
degree-and- order  12  truncation  based  on  satellite  data. 

A  perfect  spherical  harmonic  model  of  any  given  finite 
degree  and  order  is  unrealistic. 

A  study  based  on  the  thesis  of  perfect  spherical 
harmonic  modeling  is  a  valid  experiment  to  determine 
sensitivity  to  model  improvements.  On  this  premise,  a 
study  is  proposed  to  determine  the  navigation  error  due  to 
perfect  spherical  harmonic  modeling  over  a  range  of  degree- 
and-order  truncations. 

The  effects  of  this  level-of-detail  variations  in 
the  gravity  model  must  be  projected  onto  the  statistical 
model  of  the  residual  field.  A  spherical  harmonic  expan¬ 
sion  of  the  residual  field  would  have  zero  coefficients  for 
all  degrees  less  than  the  truncation  number.  With  the 
anomaly  degree  variance  model,  the  statistical  model  alter¬ 
ation  is  apparent.  This  model  was  introduced  in  Sections 
V  and  VI  and  further  discussed  in  Appendix  J.  Since  the 
anomalous  field  has  zero  coefficients  in  spherical  harmonics 
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below  the  truncation  number,  it  is  reasonable  to  model  the 
statistics  similarly. 

With  the  anomaly  degree  variance  model,  the  anomaly 
correlation  is  given  by  an  equation  from  Section  VI 
cpgg(£.£')  =  cpgg(r,r*  ,t|;) 


=  2 
n=0 


Pn(cos  4 /) 


(132) 


The  coefficient  cn  is  determined  by  an  algebraic  rule  given 
in  Appendix  J.  For  this  study, 

(  0  n  <  2  or  n  <  N 


cn  = 


A(n-l) 
(n-2) (n+B) 


n  >  N  and  n  >  2 


( 140) 


where  N  is  the  degree  and  order  of  the  spherical  harmonic 
truncation.  The  anomaly  degree  variance  model  modified 
by  (i4o)  is  the  disturbance  statistical  model  for  this 


study. 


7.2  Navigation  Error  Results 

Nested  Integrals  covariance  analyses  were  conducted 
for  truncation  levels  of  2,  3 6,  72,  and  180.  Again,  these 
studies  were  based  on  the  modified  Widnall-Grundy  naviga¬ 
tion  error  propagation  model  and  the  great  circle  trajec¬ 
tory  model.  The  truncation  level  2  is  roughly  equivalent 
to  an  ellipsoidal  model  and  represents  the  unmodified 
anomaly  degree  variance  model.  The  circular-error-probable 
results  for  all  cases  are  presented  in  Figure  18.  The 
accuracy  is  generally  better  with  each  increment.  Note, 
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Harmonic  Models  of  Limited  Degree  and  Order 
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however,  that  a  more  precise  model  does  not  guarantee 
better  system  accuracy  at  all  times.  The  dynamics  of  the 
navigation  error  propagation  model  apparently  causes  the 
ellipsoidal  model  performance  to  be  better  than  more 
detailed  models  near  the  160  minute  mission  time. 

The  accuracy  improvement  implied  by  Figure  18  compari¬ 
sons  is  desirable,*  the  cost  of  producing  the  corresponding 
models  is  not.  The  system  designer  must  weigh  both  factors 
in  deciding  on  a  level  of  model  detail.  Nested  Integrals 
supplies  a  method  to  answer  the  model  accuracy,  half  of 
the  system  design  trade-off. 
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VIII.  Nonstationary  Statistics  Demonstration 

The  statistical  models  presented  in  Sections  IV,  V, 
and  VI  are  all  based  on  an  assumption  of  homogeneity.  This 
assumption  is  instrumental  in  deriving  the  mathematical 
expressions  for  the  gravity  disturbance  stochastic  model, 
but  the  homogeneous  assumption  is  not  required  for  Nested 
Integrals  analysis.  Quite  the  contrary,  the  only  require¬ 
ment  of  Q(r.r')  is  that  it  be  integrable.  This  section  pre¬ 
sents  some  rationale  which  might  be  used  in  supporting  an 
assumption  of  nonstationary  statistics  for  the  gravity 
disturbance  field.  Then,  a  simple  nonstationary  example 
is  presented  to  demonstrate  Nested  Integrals  validity  on 
this  type  of  problem. 

First,  the  statistics  of  the  disturbance  gravity  field 
(after  the  ellipsoidal  model  has  been  removed)  vary 
regionally.  Long  has  shown  [Ref  25]  considerable  varia¬ 
bility  within  the  United  States  alone.  One  can  surmise 
from  predominant  geographic  features  like  the  mountain 
ranges,  which  run  primarily  along  north-south  lines  in  this 
country,  that  geodetic  features  are  similarly  distributed. 
The  Rice  data  [Ref  12]  from  a  survey  along  the  35- th  paral¬ 
lel  gives  some  insight  into  the  residual  field  variations 
which  might  be  anticipated  on  a  transcontinental  mission. 

The  mission  region  definition  is  the  key  to  whether 
the  nonstationary  question  needs  to  be  addressed.  If  the 
mission  definition  is  vague  and  worldwide  in  scope  (e.g. 
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for  a  short-range  transport  mission),  worldwide  average 
statistics  seem  reasonable.  If  the  mission  region  is 
restricted  (e.g.  for  MINUTEMAN  missiles),  regional  charac¬ 
teristics  should  not  be  disregarded. 

Since  Nested  Integrals  has  the  inherent  capability  to 
process  a  nonstationary  model  for  the  gravity  disturbance 
correlations,  a  simple  example  is  proposed  for  demonstra¬ 
tion  purposes.  Again,  a  proper  casting  of  the  problem  will 
yield  a  case  for  which  linear  state  space  techniques  offer 
a  crosscheck.  The  restrictions,  then,  are  to  allow  the 
verification  and  are  not  required  for  Nested  Integrals  per 
se. 

The  great  circle  trajectory  model  and  the  modified 
Widnall-Grundy  navigation  error  propagation  model  are  again 
used.  The  Table  IV  data  of  Section  IV  again  applies  with 
the  exception  of  anomaly  variance  level  which  is  explained 
below. 

8. 1  Correlation  Model 

A  nonstationary  statistical  model  gives  Q(r,r')  which 
cannot  simply  be  expressed  as  Q ( )  or  Q(tj/,a).  For  this 
example,  the  nonstationary  statistics  are  summarized  as 
anomaly  variance  which  varies  with  longitude.  One  could 
hypothesize  this  model  as 

o|(r)  *  CgU)  (141 ) 

Longitude  is  defined  by  the  design  mission  trajectory, 
hence  can  also  be  considered  a  function  of  time 
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(142) 


o|(t)  =  a|[x(t)] 

Applying  the  above  rule  to  the  linear  state  space 
gravity  disturbance  model  of  Appendix  E  will  yield  a  Q  ( t) 

O 

nonstationary  temporal  noise  strength.  For  linear  state 
space  analysis,  when  longitude  changes,  the  gravity  disturb 
ance  model  also  changes  from  one  steady-state  system  to 
another  since  P  depends  directly  on  g  .  This  quasi- 

6  B 

stationary  modeling  was  used  so  that  P  would  not  have  to 

o 

be  calculated  and  integrated. 

Applying  the  above  concept  in  a  Nested  Integrals 

analysis  requires  an  interpretation.  When  evaluating 

2 

Q(rn>r^)  does  one  use  X^  or  X^  to  calculate  g^  ?  The 
answer  is  the  same  as  for  the  question  of  which  heading 
ai  or  Si  shou-*-d  be  used  transforming  coordinates  in 

which  Q  is  expressed.  D(tn)  represents  a  driving  term  for 

• 

P  at  t=tn,  the  current  time.  The  information  entering 
through  D(  • )  at  tn  concerns  the  correlations  of  u(tn)  with 
all  past  values,  but  as  with  a,  the  information  is  paired 
( tn >  ti ) .  The  correct  interpretation  is  that  the  gravity 
correlations  at  tn  apply,  hence  jg(Xn)  should  be  used.  To 
be  rigorous,  the  order  of  arguments  in  Q  should  be  inter¬ 
preted  by  rewriting  ( 29b)  as 

D(t)  =  J\  Q[r(t)  ,r(p)  ]GT(p)  $T(t,p)  dp  (143) 

xo 

and  interpreting  (28)  modified  for  the  nonstationary  statis 


tics  of  this  example  as 

Q[>(  t)  ,r(p)  ]  =  Q[t{f,r(t)]  =  Q[t{/,X(t)] 


(144) 


(145) 


2 

The  specific  rule  used  for  o_  dependence  on  \  is 

O 

a^(X)  =  (1  +  10 X)  1800  mgal2 

D 

For  the  great  circle  trajectory,  the  anomaly  variance  begins 
at  the  same  level  as  for  the  Section  IV  study.  By  the  end 
of  the  204-minute  flight,  the  variance  level  has  increased 
by  25  percent.  This  change  was  not  selected  to  be  real¬ 
istic;  rather  the  result  should  be  distinct  from  the  Section 
IV  case  (Figure  10).  With  this  understanding,  equation  (145) 
properly  interpreted  gives  a  covariance  problem  which 
either  Nested  Integrals  or  linear  state  space  methods  can 
solve . 

8.2  Comparison  of  Results 

The  circular-error-probable  computed  by  both  Nested 
Integrals  and  linear  state  space  covariance  analyses  are 
presented  in  Figure  19  for  comparison.  The  results  again 
agree  but  not  as  well  as  the  case  of  Section  IV  which  had 
static  statistics. 

This  case  demonstrates  Nested  Integrals'  ability  to 
compute  correctly  covariance  even  in  the  case  of  non¬ 
stationary  statistics.  This  example  was  simple  and  some¬ 
what  unrealistic.  Considerable  study  would  be  required 
to  generate  a  valid  nonstationary  model.  If  the  mission 
definition  warrants  it,  and  if  empirical  data  are  available 
to  support  the  nonstationary  model,  the  Nested  Integrals 
approach  can  be  used  in  performing  the  navigation  error 
covariance  analysis. 


Tim©  (min) 


IX.  Kalman  Filter  Updates 

Any  simulation  attempts  to  describe  the  response  of 
a  system  to  some  stimuli  within  some  environment.  Nested 
Integrals,  as  a  statistical  simulation,  gives  navigation 
system  errors  due  to  gravity  disturbances  within  the  envi¬ 
ronment  of  the  trajectory  and  navigation  system.  Modern 
navigation  systems  rarely  use  inertial  instruments  alone. 
External  measurements  (e.g.  radar  position  fixes)  are  used 
to  overcome  the  low-frequency  errors  to  which  purely  iner¬ 
tial  systems  are  prone.  These  updates  are  a  part  of  the 
error  propagation  environment  which  was  not  treated  in  the 
Nested  Integrals  development  in  Section  III.  The  purpose 
of  this  section  is  to  develop  modifications  to  the  Nested 
Integrals  which  will  properly  account  for  the  effects  of 
measurement  updates  on  the  navigation  error  covariance. 

The  nature  of  the  navigation  update  must  be  specified 
in  order  to  embed  the  effects  into  the  Nested  Integrals 
algorithm.  Some  assumptions  will  be  made  on  the  types  of 
systems  and  update  methods  which  are  most  likely  to  be  used. 
First,  any  system  of  high  accuracy  will  have  linear  velocity 
or  its  equivalent  in  the  system  model.  Secondly,  the  up¬ 
date  method  will  either  be  by  Kalman  filter  methods  or  by 
closely  related  schemes  (extended  Kalman  filter  for  example) . 
A  set  of  developed  theory  exists  [Ref  15 »  325-341]  concern¬ 
ing  the  evaluation  of  filter  performance  using  a  truth  model 
to  simulate  all  those  environmental  effects  omitted  from 
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the  necessarily  limited  filter  model.  A  similar  line  of 
reasoning  will  be  followed  in  this  development. 

For  Nested  Integrals  analysis,  the  truth  model  will 
include  the  gravity  disturbance  terms.  Since  these  gravity 
error  quantities  do  not  generally  fit  the  desired  Markov 
process  mold,  no  advantage  is  gained  by  naming  them  states. 
The  G(t)u(t)  driving  term  mode  of  ( 3  )  will  be  retained. 

The  filter  is  assumed  to  contain  states  for  the  velocity 
errors,  and  the  velocity  error  states  will  normally  be 
driven  by  the  gravity  disturbance  terms.  The  geoidal 
height  N  may  be  modeled  to  drive  the  system  through  the 
system  barometric  altimeter,  if  applicable.  The  point  of 
view  taken  is  that  the  filter  model  adequately  describes 
the  propagation  of  navigation  errors.  For  an  accurate 
navigation  system,  this  assumption  is  likely  to  be  valid. 
Other  states  can  be  appended  to  the  filter  model  if  neces¬ 
sary  with  only  slight  modification  of  the  following  develop¬ 
ment. 

Once  the  algorithmic  changes  are  derived,  the  question 
of  verification  must  be  addressed.  For  this  case,  the 
linear  state  space  covariance  analysis  is  again  employed 
on  a  simple  case.  The  trajectory  and  gravity  disturbance 
models  are  selected  to  comply  with  the  linear  state  space 
limitations.  The  Nested  Integrals  results  can  then  be 
compared  to  the  linear  state  space  answers  for  verification. 
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9.1  Theory  Review 

The  objective  to  account  for  Kalman  filter  updates  can 
not  be  met  without  knowing  the  Kalman  gains.  The  update 
requires  a  system  model  and  a  measurement  model,  and  these 
models  are  used  to  calculate  Kalman  gains  in  an  analysis 
separate  and  distinct  from  the  Nested  Integrals  analysis. 
This  theory  must  be  reviewed  before  the  algorithm  changes 
are  developed. 

The  filter  model  equivalent  of  (  3 )  is 

xf(t)  =  Ff(t)xf(t)  +  Gf(t)wf(t)  (146) 

where  xf(  t)  is  an  n- vector  of  error  states  which  affect 
navigation  performance, 

Ff(t)  is  the  nXn  state  propagation  matrix, 

» 

w^Ct)  is  an  m^.- vector  of  white  Gaussian  noise,  and 
Gf(t)  is  an  nXmf  distribution  matrix. 

The  noise  w^t)  satisfies  the  following 

=  o  (147) 

<£twf(t)  w*(t+r)J  =  Qf(t)  6(r)  (148) 

The  filter  noise  model  is  not  limited  to,  indeed  may  not 
even  address,  gravity  disturbances.  The  key  point  here  is 
that  the  filter  environment  must  be  modeled  as  a  separate 
entity  from  the  Nested  Integrals  analysis.  The  model  pre¬ 
sented  in  ( 146)  ,  ( 147)  ,  and  ( 148)  is  the  standard  approach 
in  Kalman  filter  design  [Ref  15s  291-297]*  The  subscripts 
will  be  omitted  on  most  terms  from  here  forward  since  the 
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filter  model  is  assumed  to  be  sufficiently  detailed  to 
serve  as  the  Nested  Integrals  navigation  error  propagation 
model.  Exceptions  are  and  which  are  different  than 
Q  and  G  from  ( 28)  and  (  3  )  • 

For  the  filter  model  of  ( 146) ,  it  follows  that 

Pf(t)  =  F(t)Pf(t)  +  Pf(t)FT(t)  +  Gf(t)Qf(t)G^(t)  (149) 

subject  to  the  initial  condition  P^(tQ).  Also, 

^(t,^)  =  F(t)  (Mt,^)  (16a) 

Subject  to  the  initial  condition  $(t^,t^)  =  I.  These  equa¬ 
tions  are  fundamental  in  describing  the  behavior  of  the 
filter  error- covariance  estimate  between  measurement  up¬ 
dates.  The  effects  of  measurement  updates  are  assumed  to 
be  discrete  corrections  to  navigation  estimates,  and  the 
discrete  update  of  the  Pf~matrix  is  needed  to  describe  the 
filter  operations.  Continuous  feedback  of  revised  esti¬ 
mates  can  also  be  treated  with  a  few  modifications  to  the 
development  here  [Ref  150333' 

The  Nested  Integrals  solution  is  based  on  a  discrete 
solution  time  sequence  tQ  <  t^  <  • » •  <  t^,  symbolically 
(tn).  The  measurement  updates  are  assumed  to  come  at  time 
points  which  are  a  subsequence  of  this  time  base.  That  is, 

measurements  are  at  t  <  t  <  • • •  <  t  ,  symbolically 

n  1  n2  M 

(t  ).  At  each  time  t  ,  an  r-dimensional  vector  measure- 
ni  ni 

ment  occurs  and  is  described  by 


z(tni)  =  H(tni)x(tn_)  +  v(tn:) 


(150) 
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where 


H(tn  )  is  an  rXn  measurement  matrix,  and 

v(t  )  is  an  r-dimensional  discrete,  white,  Gaussian 


noise  process  satisfying 


<?  [-(V]  =  0 

^[Y(tn.)yT(tn3)J 


R(t  )  i=J 
ni 


(151) 


(152) 


Let  t  indicate  the  condition  of  a  quantity  prior  to 
ni 

the  update  at  t-  ,  and  let  t*  represent  the  post-update 
ni  ni 

conditions.  Then,  using  (149), 

pf(t:  )  =  $(t;  ,  t+  )  ?  (t+  )  $T(t;  ,t*  ) 

i  ni  ni  ni_i  1  ni_i  ni  ni_i 

+  J  1  . •P)Gf(p)Qf(p)Gj(p)$T(t"  iPldP 

^i.l  1  1  (153) 

at  tn  the  Kalman  gain  is  calculated  by 

k<\.>  ■  >  r«<vpf(vHT<v +  R(v,r 

i  i  ii_i  i  i  i  j 

(154) 


The  update  of  the  filter  covariance  matrix  is  given  by 
>  =  pf(^  >  -  K(tn  >»<*„  > 


l  "i 


♦  K(t  )R(t  )KT(t  ) 
1  "i  ^ 


(155) 
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The  equation  set(l6a),  (153.) »  (1 5*0  »  and  (155)  allow  the 
computation  of  the  Kalman  filter  gains  K(tn  )  needed  to 
update  the  covariance  calculated,  separately,  by  the  Nested 
Integrals  method. 


9.2  Nested  Integrals  Algorithm  Changes 

The  Nested  Integrals  algorithm  calculates  the  covari¬ 
ance  between  updates  correctly.  That  is,  given  P(t  ) 

ni-l 

and  for  t  <  t  <  t ,  equation  (29a)  is  valid 
ni-l  ni 

P(t)  =  F(t)P(t)+P(t)FT(t)+G(t)D(t)+DT(t)GT(t)  (29a) 


But,  since  D(t)  inherently  includes  the  propagation  of 
errors  in  its  definition,  the  effects  of  updates  on  D(t) 
must  also  be  included  in  the  Nested  Integrals  algorithm. 
Recall 


D( t)  =  J\  Q(t,p)GT(p)  $>T(t,p)dp 


(29b) 


To  reflect  the  measurement  update  on  both  P(t)  and  D( t) , 
the  error  state  x  is  propagated  across  the  update  as 

x(t+  )  =  [I-K(t  )H(t  )3x(t:  )  +  K(t  )v(t  )  (156) 
ni  ni  ni  ni  ni  nA 

But  v(t  )  is  zero  in  this  gravity-error-only  analysis,  and, 
ni 

this  form  can  be  used  to  interpret  a  discrete  state  propa¬ 
gation  matrix  as 

$>(<  ,t"  )  =  I-K(  t  )H(  t  )  (15?) 

ni  ni  ni  ni 
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Since  Nested  Integrals  involves  geodetic  errors  alone,  the 
R-matrix  of  (155)  is  not  involved  in  the  P(t)  update.  One 
should  note  that  any  altimeter  errors  were  assumed  to  enter 
the  error  process  through  a  u(t)  element.  One  could  pro¬ 
pose  altimeter  updates  in  which  case  v  would  be  correlated 
with  U{  this  approach  is  not  taken  here.  Using  equation 
(155)  as  a  guide,  the  covariance  update  is 


p<< .)  '  >p(tn.)$T<tn.-tn  > 

1  111  11 


(158) 


The  update  of  D(t)  must  be  derived  since  an  equivalent 
matrix  is  not  normally  treated  in  Kalman  filter  theory. 

The  integrand  of  (29b)  is  assumed  to  be  finite,  so  the 
upper  limit  of  the  integration  can  be  moved  from  one  side 
of  the  update  to  the  other  without  destroying  the  equality. 
That  is,  the  interval  (t_,t+)  is  assumed  to  be  a  set  of 
zero  measure.  The  resulting  D(t)  across  updates  can  be 


found  by 


9<tV 


p)GT(p)  $T( t*. ,p)dp 


{t“ 

1  Q( t"  ,p)GT(p)  $T(t*  , p) dp 

t  ni  ni 

o 

<t“ 

t  1  Q^n.  ^(t^  ,t”  )$(t"  ,p)J  Tdp 


1  Q(t->lP)GT(P)  $T(t"  ,p)dp  UT(t+  ,t"  ) 

m  1  1  J  ni  ni 


=  D(t- )  ^T(t:  ,t- ) 


ni  ni 


(159) 


This  update  for  D(t)  will  force  discontinuities  at  t  . 

ni 

In  solving  the  differential  equation  ( 29a)  for  the  inter¬ 
vals  between  updates,  the  pre-  or  post-update  D(t)  must  be 
selected  with  care.  Let  j  =  n. j  for  the  interval  t.  , 
to  t-,  the  P  calculation  should  be  based  on  D(t.  ,)  and 

t)  V  "L 

D(tT).  When  t.  is  not  a  measurement  update  time,  the  t. 

J  J  J 

and  t~  data  are  the  same . 

J 

The  discrete  error  propagation  defined  by  (  15?) ,  there¬ 
fore,  gives  the  key  to  the  Nested  Integrals  Kalman  filter 
update.  Using  the  result  of  (1 57),  P(t)  is  updated  as 
specified  in  (158)  and  D(t)  is  updated  as  specified  in 
(159)-  The  only  new  information  required  is  the  Kalman 

gain  K(t  )  and  the  measurement  matrix  H(t  ). 
ni  ni 


9.3  Verification  Case 

Linear  state  space  techniques  are  the  natural  mode  of 
Kalman  filter  design  and  analysis.  A  linear  state  space 
covariance  analysis  can  be  defined  to  include  Kalman  filter 
updates,  and  the  results  used  as  a  benchmark  to  judge  the 
modified  Nested  Integrals  algorithm.  The  great  circle  case 
of  Section  IV  will  again  serve  as  the  vehicle  for  comparing 
the  results  from  the  two  covariance  analysis  methods.  The 
measurement  model  and  the  filter  model  must  be  specified  in 
addition  to  the  three  models  (trajectory,  error  propagation, 
and  disturbance  correlation)  of  the  great  circle  case. 

Some  additional  explanation  is  required  concerning  how  the 
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linear  state  space  method  is  modified  to  incorporate  the 
effects  of  the  Kalman  filter  update. 


The  two  new  models  required  are  described  below.  The 
linear  state  space  covariance  method  is  described  in 
Section  II;  recall  that  the  covariance  of  an  augmented  state 
is  calculated.  The  augmented  state  is  given  by 
x(t) 


Vt)! 


(7  ) 


The  effect  of  the  update  equation  above  is  on  the  x(t) 
states  alone.  Using  the  notation  for  discrete  error  propa¬ 
gation  matrices, 

0 


I-K(t  )H(t  )  ; 

1  1  I 

_ J _ 

o  !  o 


(160) 


The  linear  state  space  augmented  state  covariance  update 
is  given  by 


.t-  )pQ(t:  )$?(<  ,t;  ) 


a'  n^ 


a'  n^'  ni'  a'  n^'  a'  n^’  ni 


(  161) 


For  the  time  intervals  between  updates,  equation  (29)  will 
describe  the  dynamics  of  P(t).  The  modification  of  the 
Section  II  described  linear  state  space  covariance  analysis 
by  (161)  completely  specifies  the  analysis  by  that  method. 

A  practical  note  should  be  made  before  proceeding  with 
the  model  definitions.  The  predictor-corrector  integration 
algorithm  used  to  solve  (14)  and  (29a)  is  not  designed  to 
handle  discontinuities.  The  discontinuities  at  the  update 
times  would  create  large  predictor  errors  for  which  the 
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computed  corrections  might  he  erroneous.  This  problem  is 
avoided  if  one  reinitializes  the  integration  algorithm 
after  each  update.  The  instructions  for  this  initializa¬ 
tion  are  a  normal  part  of  the  software  description  with 
such  algorithms. 


9.3.1  Filter  Model 

To  keep  the  verification  case  relatively  simple,  no 
additional  states  are  considered  other  than  the  nine  states 
of  the  modified  Widnall-Grundy  navigation  error  model. 
Examples  of  possible  additional  states  are  gyroscopic  drift 
and  accelerometer  bias.  The  system  noise  model  simulates 
the  effects  of  such  errors  through  white  noise  driving  the 
velocity  (x^,x^,  and  x^)  and  attitude  (x^,Xg,  and  x^) 
derivatives.  Interpreting  this  filter  model  choice  in  terms 
of  equation  (146),  define 

Wf(t)  =  <  :  >  (162) 

[w6(t)J 

and  the  associated  distribution  matrix  is 


Gf(t)  =  Gf  = 


(163) 


where  G^.  is  9X6  and  the  identity  matrix  is  6X6.  The  noise 
strength  matrix  associated  with  w^.  is  given  by 
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(164) 


The  acceleration-level  white  noise  is  specified  at  a  level 
of  Qy  =  16x10”^  ft^/sec-^,  and  the  attitude-rate-level  noise 
is  modeled  as  =  36x10 ~  radian  /sec.  These  noise 
levels  define  an  inertial  navigation  system  of  approximately 
one  nautical  mile  per  hour  circular-error-probable  growth 
rate . 


The  modified  Widnall- Grundy  x('t)  ant*  F(t)  along  with 
(162)  through  (164)  give  the  filter  structure.  This  struc¬ 
ture  driven  by  the  noise  model  above  governs  the  filter 
model  covariance  simulation. 


9-3*2  Measurement  Model 

A  horizontal  position  fix  is  hypothesized  as  the 
measurement  external  to  the  inertial  navigation  system. 

A  radar  system  with  a  100-foot  root-mean-square  resolution 
takes  east  and  north  position  coordinate  measurements  at 
intervals  throughout  the  mission.  For  simplicity,  the 
measurements  are  made  on  five-minute  intervals.  This  rate 
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corresponds  to  once  every  ten  Nested  Integrals  integra¬ 
tion  steps  —  recall  the  30-second  step  size  from  Table  IV. 

The  horizontal  position  measurements  allow  the  obser¬ 
vation  of  east  position  error  xg  and  north  position  error 
as  described  belowj 


xg  =  r  cos0  6X  =  r  cos0  x-^ 
xn  =  r  60  =  r  x2 
With  (165)  and  (166) ,  define 


and  for  (150) 


(165) 

(166) 


(167) 


tr  cos0  00000000 
0  rOOOOOOO 

The  noise  strength  of  (152)  is  modeled  as 


(168) 


R  = 

2  2 

where  or  is  (100  ft)  by  the  above  assumption. 

The  measurement  model  is  complete  with  ( 168)  and  ( 169) 
providing  the  necessary  components  of  ( 154)  and  ( 155)  . 


(169) 


9.4  Comparison  of  Results 

Based  on  the  filter  and  measurement  models,  a  set  of 
Kalman  gains  were  calculated  and  stored  for  use  in  both  the 
linear  state  space  and  Nested  Integrals  covariance  analyses. 
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Then,  separate  results  were  generated  using  these  alterna¬ 
tive  statistical  methods.  The  circular-error-probable 
for  both  methods  is  plotted  in  Figure  20.  The  two  curves 
are  practically  indistinguishable,  and  the  curve  markers 
are  again  offset  to  clarify  that  two  curves  are  plotted. 

The  conclusion  from  this  comparison  is  that  the  modi¬ 
fications  to  Nested  Integrals  can  be  used  to  account 
properly  for  the  effects  of  Kalman  filter  updates  through¬ 
out  a  mission.  As  in  Section  IV,  the  linear  state  space 
solution  is  known  to  be  correct  on  this  restricted  trajec¬ 
tory.  Nested  Integrals  has  the  capability  to  compute  co- 
variance  on  a  wider  variety  of  missions  and  with  more 
flexibility  in  the  choice  of  the  gravity  disturbance 
statistical  model,  and  now  has  been  extended  to  include  the 
flexibility  of  employing  Kalman  filter  updates  on  the 
navigation  mission. 
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The  objective  of  this  research  was  to  develop  a  new 
computational  technique  which  would  provide  a  performance 
assessment  of  alternative  gravity  models  on  realistic 
scenarios.  This  section  provides  a  summary  of  the  develop¬ 
ment  of  the  Nested  Integrals  algorithm.  The  verification 
cases  and  capability  demonstrations  are,  also,  recapitu¬ 
lated  for  completeness. 

A  navigation  system  accuracy  measure  was  selected  as 
the  model  comparison  performance  measure.  Statistical 
methods  of  producing  such  a  measure  were  selected  since 
these  methods  admit  the  limitations  in  modeling  the  gravity 
field  and  permit  analysis  on  a  general  mission  definition. 
Two  present  statistical  methods  were  considered  and  dis¬ 
cussed-.  linear  state  space  covariance  analysis  and  Monte 
Carlo.  The  covariance  integral  was  developed  as  a  com¬ 
promise  offering  more  flexibility  than  linear  state  space 
covariance  analysis  but  promising  a  savings  in  computa¬ 
tional  costs  compared  with  Monte  Carlo. 

Numerical  methods  were  developed  to  compute  navigation 
error  covariance  based  on  an  integral  expression  which  re¬ 
lates  the  covariance  to  the  gravity  disturbance  causes. 
Three  models  are  required  for  these  analyses.  A  linear, 
first-order  differential  equation,  ( 3) »  which  relates 
the  navigation  error  rates  to  navigation  errors  and  to  the 
gravity  disturbances.  Because  inertial  error  propagation 
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is  trajectory  dependent  and  because  the  gravity  disturbance 
is  a  spatial  function,  a  design-mission  trajectory  model  is 
needed  to  represent  the  dynamics  and  geometry  of  an  ensembl 
of  missions  which  form  the  statistical  sample  base. 

Finally,  the  statistics  of  gravity  disturbances  encountered 
over  this  ensemble  of  missions  must  be  summarized  in  the 
form  of  the  correlation  model.  This  correlation  function 
is  the  key  element  in  the  integrand  of  the  covariance 
integral. 

Three  separate  numerical  algorithms  were  developed 
to  approximate  the  covariance  integral.  Two  of  these  were 
direct  integral  approximations;  one  used  a  rectangular 
rule  and  the  other  trapezoidal  in  approximating  the  inte¬ 
gral.  The  third  algorithm  was  created  by  taking  the  time 
derivative  of  the  covariance  integral  using  Leibnitz'  rule, 
equation  (29).  The  covariance  integral  was,  thereby, 
recast  from  a  double  integral  into  two  single  integrals 
which  are  nested.  The  outer  of  these  Nested  Integrals  is 
in  a  form  compatible  with  standard  predictor-corrector 
numerical  solution  techniques. 

To  implement  any  of  these  three  numerical  methods 
would  be  computationally  very  costly  if  a  naive  approach 
were  taken.  To  avoid  the  requirement  to  store  state 
transition  matrices  between  each  pair  of  evaluation  points, 
the  semi-group  property  of  state  transition  matrices  was 
used.  This  reduced  the  data  storage  requirement  to  approxi 
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mately  the  square  root  of  the  storage  required  by  the 
straightforward  approach.  To  accommodate  this  storage 
savings  and  avoid  data  retrieval  costs,  an  indexed  sequen¬ 
tial  file  structure  was  selected  for  the  required  state 
transition  matrix  and  position  data.  This  approach  permits 
relatively  efficient  sequential  file  use  once  the  initial 
record  is  located.  To  exploit  this  feature,  summations 
which  occur  in  the  integral  approximations  are  rim  in  the 
reverse  chronological  order.  Together,  the  file  structure 
and  algorithmic  form  give  an  economically  viable  method  of 
computing  covariances. 

Of  these  three  methods,  Nested  Integrals  proved  to 
be  more  efficient  and  of  equal  accuracy.  The  analysis  of 
an  undamped  Schuler  loop  driven  by  exponentially  correlated 
noise  was  used  as  a  medium  for  this  tradeoff.  This  example 
problem  demonstrated  graphically  the  need  to  select  inte¬ 
gration  step  size  small  enough  to  meet  the  Shannon  rate 
and  correctly  represent  the  frequency  content  of  the  inte¬ 
grand  of  the  covariance  integral.  Results  on  the  subse¬ 
quent  damped  Schuler  loop  example  demonstrated  that  the 
error  level  is  lower  in  comparison  since  the  system  model 
is  inherently  stable. 

Once  the  Nested  Integrals  method  was  selected,  a  full- 
scale  verification  study  was  conducted  using  a  great  circle 
trajectory  and  a  shaping  filter  gravity  disturbance  model. 
The  verification  was  successful,  and  the  same  comparison 
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format  was  used  to  demonstrate  the  errors  in  linear  state 
space  covariance  analysis  on  more  complex  trajectories. 
Minor  circle  trajectories  were  used  to  demonstrate  that 
linear  state  space  results,  while  correct  on  a  great  circle 
case,  are  increasingly  in  error  as  the  circle  radius  is 
decreased. 

Next,  the  Nested  Integrals  method  was  verified  against 
an  independent  Monte  Carlo  study  of  a  highly- dynamic  air- 
launched  missile.  This  study  represents  a  real-world 
application  of  Nested  Integrals.  The  Monte  Carlo  study 
models  were  emulated  in  the  Nested  Integrals  analysis.  For 
this  study,  a  different  functional  form  for  the  correla¬ 
tion  model  was  specified  in  the  Monte  Carlo  analysis,  and 
the  statistical  model  flexibility  inherent  in  the  Nested 
Integrals  formulation  was  required.  The  results  from  both 
methods  compared  well  and  offer  an  additional  validation 
of  Nested  Integrals.  For  this  problem,  Nested  Integrals 
demonstrated  a  30:1  computational  time  advantage  over 
Monte  Carlo.  This  efficiency  was  a  prime  objective  of  the 
Nested  Integrals  method. 

Clearly,  Nested  Integrals  meets  the  original  expec¬ 
tations:  it  is  a  more  flexible  analysis  technique  than 

linear  state  space  covariance  analysis,  and  it  is  computa¬ 
tionally  more  efficient  than  Monte  Carlo.  Three  demonstra¬ 
tions  were  performed  to  give  further  examples  of  specific 
Nested  Integrals  capabilities.  Three  different  correlation 


models  were  analyzed  using  the  same  trajectory  and  error 
propagation  model  for  each  study.  This  study  demonstrated 
the  correlation  model  flexibility  of  Nested  Integrals  and, 
simultaneously,  demonstrated  the  variability  in  results 
possible  with  existing,  supposedly  equivalent,  correlation 
models..  Next,  the  ability  to  compare  gravity  models  was 
demonstrated  using  perfect  spherical  harmonic  modeling  as 
a  case  study.  The  final  demonstration  verified  the  capa¬ 
bility  of  Nested  Integrals  to  account  properly  for  non¬ 
stationary  statistics,  a  correlation  model  nuance  which 
might  have  application  on  some  missions. 

In  a  final  development,  modifications  to  Nested 
Integrals  algorithm  were  made  to  accommodate  Kalman  filter 
updates.  This  capability  is  needed  to  analyze  most  modem 
navigation  missions.  A  linear  state  space  covariance 
analysis  on  the  great  circle  trajectory  was  used  to  verify 
the  modified  Nested  Integrals  method. 

With  these  developments,  demonstrations  and  verifi¬ 
cations,  a  new  statistical  analysis  technique  has  been 
presented.  This  technique  offers  a  more  flexible  alterna¬ 
tive  to  the  common  linear  state  space  method  while  being 
more  computationally  efficient  than  Monte  Carlo.  The 
Nested  Integrals  algorithm  is  a  general  method  for  producing 
the  navigation  system  accuracy  statistical  performance 
index  by  which  gravity  models  can  be  compared. 
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XI .  Recommendations  for  Future  Research 

This  research  has  produced  a  gravity  model  evaluation 
method  which  can  be  applied  on  general  mission  scenarios. 

This  section  contains  some  suggestions  to  improve  this 
statistical  method  and  suggestions  of  possible  implemen¬ 
tations  which  require  additional  research. 

First,  to  develop  this  method  further,  it  is  suggested 
that  the  Nested  Integrals  algorithm  be  developed  for  auto¬ 
mated  variable  step  size  selection.  The  trajectory  model 
would  have  to  interact  directly  with  Nested  Integrals. 

Since  approximation  error  is  usually  proportional  to  a 
power  of  integration  step  size  on  simple,  one-dimensional 
integration  rules,  one  might  calculate  the  coefficient  at 
each  point  and  select  the  next  step  size  based  on  the 
reciprocal  of  this  norm.  This  time-step  calculation  would 
tend  to  give  constant  potential  for  approximation  error. 

Large  time  steps  would  be  allowed  when  the  coefficient  is 
small,  and  during  highly  dynamic  regions  the  coefficient 
would  be  large  forcing  smaller  time  steps.  The  end  result 
could  be  a  substantial  increase  in  efficiency  since  compu¬ 
tational  costs  are  roughly  quadratic  with  number  of  points. 

A  comprehensive  study  of  the  accuracy  bounds  of  the 
Nested  Integrals  algorithm  would  be  useful  in  future  appli¬ 
cations.  One  approach  would  be  to  treat  this  as  a  frequency 
response  problem.  The  correlation  function  could  be  cos  (oijf 
where  co  is  the  frequency  parameter.  Using  this  correlation 
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analyses  can  be  conducted  over  a  range  of  a)  values.  The 
trajectory  and  error  propagation  models  must  be  given.  A 
single  error  quantity,  like  maximum  circular-error- 
probable,  plotted  versus  (i)  would  be  the  form  of  output. 
Another  parameter  in  this  study  would  be  integration  step 
size.  These  results  would  be  trajectory  dependent,  and  the 
true  covariance  would  be  required  to  make  algorithm  error 
observable. 

The  batch-processing  nature  of  Nested  Integrals  will 
limit  its  use  in  real-time  applications.  As  a  post¬ 
mission  analysis  tool,  it  offers  some  attractive  features. 

For  test  flights  of  navigation  systems,  Nested  Integrals 
can  be  used  to  account  statistically  for  gravity  induced 
errors  on  complex  trajectories.  Better  estimates  of  instru¬ 
ment  errors  would  be  the  motivation  for  this  analysis. 
Research  is  required  to  incorporate  Nested  Integrals  into 
the  post  flight  smoothing  and  error  identification  processes. 

Along  a  similar  line,  Nested  Integrals  should  prove 
useful  in  gravity  surveys  conducted  with  inertial  naviga¬ 
tion  for  position  reference.  The  algorithm  should  prove 
useful  in  post-survey  data  processing.  With  adaptation, 
the  method  developed  here  should  also  be  a  useful  tool  in 
survey  design. 

Another  application  is  to  evaluate  different  gravity 
modeling  techniques.  The  example  presented  in  Section  VII 
is  simple;  any  realistic  comparison  would  pose  a  distinct 
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problem  for  the  analyst.  How  does  one  account  for  gravity 
model  effects  on  the  residual  field?  Models  can  be  of 
different  functional  form  (e.g.  global  spherical  harmonics 
versus  local  Chebychev  polynomial  functional  fits)  and 
of  different  levels  of  detail  within  one  model  type  (e.g. 
different  order  Chebychev  polynomials  within  the  same 
region) .  Finding  a  general  framework  within  which  both 
types  of  problems  can  be  solved  is  a  substantial  research 
task.  The  residual  field  is  treated  as  an  error  field  with 
first-order  approximations  throughout  the  theory.  Perhaps 
the  gravity  model  improvements  can  be  represented,  to 
first-order,  as  a  bounded  linear  operator  over  this  error 
field.  Symbolically 

T'(r)  =  {?(*)} 

where  T'  is  the  residual  potential  field  after  the  improved 
gravity  model  and  operates  on  T,  the  residual  field 
from  the  ellipsoidal  model.  The  covariance  of  the  residual 
field  would  be  given  by 

> 

where  the  operator  subscripts  indicate  which  position  vector 
would  be  involved.  This  idea  would  require  a  substantial 
research  effort  to  bring  it  to  fruition. 

Research  is  needed  on  two  issues  concerning  the 
statistical  model  for  gravity  disturbance  correlations. 
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First,  although  all  three  of  the  modeling  methods  studied 
in  Section  VI  are  intended  to  approximate  the  empirical 
correlation  function,  there  is  no  mention  of  completeness 
of  the  underlying  approximations.  That  is,  "can  any  of 
the  available  modeling  options  replicate  the  empirical 
correlation  function  within  some  arbitrarily  small  norm?" 
Second,  the  analyses  conducted  herein  assumed  the  residual 
field  statistics  can  be  "upward  continued"  [[Appendix  A 
and  Ref  13]  from  some  reference  surface.  What  if  the 
inertial  navigation  system  (equivalent  if  not  explicit  ) 
gravitation  model  is  not  harmonic?  Such  approximations 
exist  in  less  precise  inertial  systems  than  addressed  herein. 
With  these  existing  non-harmonic  models,  how  can  correla¬ 
tions  of  gravity  errors  aloft  be  surmised  from  reference 
surface  statistics?  For  both  of  these  potential  research 
areas,  the  first  question  to  be  answered  is  whether  or  not 
these  statistical  modeling  questions  are  important  in  terms 
of  the  potential  impact  on  answers.  That  is,  what  is  the 
sensitivity  of  navigation  error  statistics  to  empirical 
correlation  function  approximation  errors  or  to  non¬ 
harmonic  gravitation  model  effects. 

Kalman  filter  effects  were  included  in  the  Nested 
Integrals  algorithm  in  Section  IX.  The  measurement  update 
tends  to  sever  the  correlation  of  present  navigation 
errors  to  past  gravity  field  errors.  The  effects  of  mis- 
modeling  the  long- time-range  gravity  disturbance 
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correlations  in  Q(r,r')  should  be  less  important  because 
of  the  update.  A  navigation  system  with  frequent  updates 
might  be  adequately  treated  by  the  more  efficient  linear 
state  space  method  even  with  trajectory  restriction  viola 
tions.  Numerical  studies,  such  as  the  minor  circle  para¬ 
metric  study  of  Section  IV,  could  give  insight  into  this 
issue.  A  comprehensive  study  would  be  useful  in  planning 
the  statistical  analysis  of  future  systems  and  missions. 

These  research  areas  seem  the  most  promising  exten¬ 
sions  and  applications  of  Nested  Integrals.  Any  number 
of  small  studies  similar  to  Sections  VI,  VII,  and  VIII 
of  this  work  would  contribute  to  understanding  the  uses 
and  limits  of  Nested  Integrals.  Additional  comparisons 
between  Nested  Integrals  and  linear  state  space  analysis 
would  be  useful  in  further  defining  the  accuracy  limits 
of  the  efficient  linear  state  space  method. 


APPENDIX  A 


Comments  on  Gravity  Correlation  Models 

A  gravity  disturbance  correlation  model  is  required 
in  order  to  perform  the  analyses  which  are  the  subject  of 
this  research.  This  model  conveys,  in  an  average  sense, 
the  interrelationships  of  all  gravity  disturbance  terms 
which  affect  the  navigation  accuracy.  The  choice  of  trajec¬ 
tory  and  navigation  error  propagation  models  is  straight¬ 
forward  in  comparison  to  the  choice  of  this  statistical 
model.  No  panacea  exists  in  this  case.  A  universally 
accepted  model  which  can  be  used  in  all  cases  has  not  been, 
and  probably  cannot  be,  defined.  The  cases  which  might 
create  problems  for  a  general  model  include  geographically 
restricted  missions  and  gravity  model  improvements.  In 
either  of  these  examples,  a  correlation  model  must  be  formed 
in  a  manner  peculiar  to  the  problem  being  addressed.  The 
purpose  of  this  appendix  is  to  illuminate  some  of  the 
issues  which  should  be  considered  in  choosing  the  gravity 
disturbance  correlation  model. 

The  first  problem  an  analyst  faces  when  surveying  model 
concepts  is  language.  The  fact  that  few  orderly  threads 
weave  through  this  labyrinth  of  model  types  compounds  the 
confusion.  No  clear  taxonomy  classifies  these  models  for 
the  novice,  so  a  few  of  the  labels  applied  to  models  are 
defined  as  follows: 
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a.  Homogeneity.  [Ref  4:853  The  statistics  are  not 
a  function  of  position  (i.e.  spatially  stationary 
statistics) .  Obviously,  gravity  quantities  change 
with  altitude,  so  this  term  describes  the  behavior 
at  one  altitude  level. 

b.  Isotropicity .  [Ref  4s 853  The  statistics  are 
independent  of  direction  or  heading  in  the  hori¬ 
zontal  plane. 

c.  Global.  The  statistics  are  worldwide  averages 
or  expectations  as  opposed  to  local  or  regional. 

d.  Anomaly.*  The  statistics  are  based  on  gravity 
anomaly  Ag-  The  symbol  cp  ( ■ )  is  used  in  this  work 

6o 

for  anomaly  correlation.  Moritz  [Ref  4:833  and  others 
in  Geodesy  use  C(-)  for  this  function. 

e.  Potential.*  The  statistics  are  of  or  based  on 
anomalous  potential  T.  The  symbol  cpTT(  • )  is  used 
herein;  others  [Ref  4:863  use  K( • )  for  this  function. 

f.  Undulation.*  The  statistics  are  of  or  based  on 
undulation  N  of  the  geoid  (also  known  as  geoidal 
height) .  This  quantity  is  directly  related  to 
anomalous  potential  T  by  Brun’s  formula  [Ref  33 

T(r)  =  g  N(r)  (A-l) 

Therefore,  the  correlations  are  related  by 

CpTT(  ‘ )  =  g2  • )  (A-2) 

*  Definitions  for  Ag,  N,  T,  and  other  gravity  disturbance 
quantities  are  given  in  Reference  3. 
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g.  k-th  order.  This  term  states  the  order  of  Markov 
process  used  in  the  model.  For  example,  a  third- 
order  undulation  model  is  based  on  undulation 
correlation  and  structured  as  a  third-order  Markov 
process. 

h.  Upward  Continuation.  Not  a  classification  per  se, 
this  term  refers  to  the  extrapolation  of  statistics 
aloft  from  the  datum  surface.  Some  models  (e.g.  the 
Tscherning-Rapp  anomaly  degree  variance  model  [Ref 
2l3  in  Sections  V  and  VI)  have  this  inherent  capa¬ 
bility;  other  models  do  not  (e.g.  the  linear  state 
space  model  [Ref  133  introduced  in  Section  IV) . 

i.  Self-consistent.  The  model  generates  gravity 
disturbance  term  statistics  (e.g.  auto-  and  cross¬ 
correlations)  which  are  consistent  with  respect  to 
the  gravitational  field  theory  from  which  interrela¬ 
tionships  of  these  terms  can  be  derived. 

These  terms  and  more  await  the  researcher  delving  into  this 
modeling  issue.  The  three  full-scale  models  discussed  in 
Section  VI  [Refs  19,  20,  21,  and  22]  offer  a  starting  point 
at  least.  References  2,  4,  and  12  offer  some  tutorial 
assistance,  also.  In  the  end,  whether  or  not  a  model  is 
adequate  will  be  the  subjective  judgment  of  the  analyst, 
made  either  consciously  or  by  default. 

Consider  now  a  scenario  of  how  a  model  might  be  formed. 
Empirical  correlation  functions  can  be  generated  from  gravity 
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Figure  A-l.  Development  of  a  Gravity  Disturbance 
Statistical  Model 
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(or  gravity  related)  measurements.  A  model  structure  is 
next  either  selected  from  existing  examples  or  derived  on 
some  new  rationale.  The  model  parameters  are  selected  to 
replicate  the  empirical  function.  And  finally,  the  mathe¬ 
matical  structure  and  the  data  parameters  of  the  model  are 
programmed  as  a  subroutine  which  accepts  position  pairs  as 
input  and  which  yields  the  required  correlations  as  outputs. 
The  correlations  or  at  least  representatives  should  be 
compared  to  the  empirical  function  for  validation. 

The  first  step  assumes  gravity  data  is  available.  The 
primary  data  available  for  such  analysis  are  gravity  magni¬ 
tude  measurements.  These  data  are  converted  to  free  air 
anomaly  and  referenced  to  the  geoid.  So,  any  averages 
formed  from  these  data  apply  to  the  geoid  or  an  approxi¬ 
mating,  analytical  surface  like  the  Bjerhammer  sphere  CRef 
3  i32l].  These  geodetic  data  processing  steps  introduce 
error  through  the  models  CRef  33  used  to  produce  the  anomaly 
and  to  estimate  the  value  on  the  geoid.  The  data  produced 
by  such  processing  is  no  longer  truly  "empirical".  Corre¬ 
lation  models  based  on  these  processed  data  will  be  referred 
to  as  "empirical",  but  the  reader  should  be  aware  of  these 
possibly  corrupting  influences. 

Deflections  of  the  vertical  have  been  measured  by 
astro-geodetic  means  and  geoidal  undulations  by  satellite 
altimetry.  These  data  could  also  be  processed  to  form 
empirical  correlation  functions.  Indeed,  heterogeneous 


types  of  data  can  be  used,  if  coverage  is  adequate,  to 
produce  cross-correlations  as  well.  Such  steps  need  not 
be  taken  for,  as  will  be  shown  later,  only  one  correlation 
function  is  needed.  The  other  auto-  and  cross-correlations 
can  be  derived  from  this  base  function  using  interrelation¬ 
ships  from  gravitational  theory  which  is  assumed  to  correctly 
represent  the  physics  of  this  problem. 

Data  for  model  development  is  limited  by  economics 
and  by  politics;  therefore,  compromises  must  be  made  in 
representing  the  correlation  function.  The  manner  in  which 
the  empirical  correlation  is  formed  makes  a  statement  about 
the  type  of  model  it  will  support.  If  the  correlation  is 
calculated  by  averaging  data  over  the  earth's  surface,  the 
resulting  correlation  function  will  represent  global  sta¬ 
tistics.  If  the  independent  variables  for  the  averaging 
are  relative  position  quantities,  the  average  represents 
a  model  based  on  the  homogeneous  assumption  [Ref  4  :85]. 

If  the  independent  variable  is  central  angle  (or  equivalent 
surface  arc  distance)  shift  with  the  average  being  over  all 
headings,  the  resultant  empirical  correlation  is  consis¬ 
tent  with  the  isotropic  assumption  [Ref  4  s  85 U •  With  all 
of  this  structure,  the  reality  of  the  data  availability 
is  that  global  high  frequency  coverage  does  not  exist  due 
to  the  magnitude  of  the  task  implied.  The  homogeneous 
and  isotropic  assumptions  will  undoubtedly  be  made  just 
to  bring  more  data  points  into  the  empirical  formulation. 

1?8 


Assuming  that  the  empirical  correlation  exists,  the 
model  structure  must  be  defined.  The  examples  in  Section 
VI  are  fully-developed  and  appear  adaptable.  Whether  these 
models  can  fit  an  arbitrary  empirical  correlation  function 
with  arbitrarily  small  error  is  a  question  worthy  of  atten¬ 
tion.  For  now,  assume  the  model  can  be  fit  to  the  empirical 
function  with  acceptable  veracity.  Then,  one  merely  trans¬ 
forms  the  mathematical  model  into  a  computer  subroutine 
and  proceeds  with  the  analysis. 

The  question  may  be  raised,  "Why  go  to  these  mathemati¬ 
cal  forms  when  the  empirical  function  is  the  desired  corre¬ 
lation?"  Obviously,  a  table  look-up  approach  can  be  used 
on  the  basic  correlation  function.  The  Q(r,r')  matrix  func¬ 
tion  of  (28)  requires  a  set  of  auto-  and  cross-correlation 
functions,  however,  and  this  set  should  be  consistent.  Data 
does  not  exist  to  produce  adequate  empirical  models  f c r  all 
of  these  gravity  disturbance  terms,  so  producing  them  indi¬ 
rectly  from  a  basic  model  is  required.  The  number  of  Q(r,r’) 
evaluations  is  great  on  even  a  modest  analysis.  The  cases 
in  Section  IV  were  based  on  409-point  Nested  Integral  analy¬ 
sis,  and  each  case  required  83,845  Q(r,r')  evaluations,  for 
instance.  A  closed  mathematical  expression  is  a  practical 
necessity  in  such  analyses. 

Fortunately,  a  closed  mathematical  expression  is 
only  required  for  the  basic  correlation  function.  The 


other  correlations  can  be  derived  from  this  source.  The 
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anomaly  or  the  anomalous  potential  correlation  is  suffi¬ 
cient  [Ref  4:94-983  as  a  base  fox  this  derivation.  To 
see  this  fact,  consider  the  following. 

The  gravity  disturbance  vector  is  related  to  potential 

by 

6g<r)  =  ^  T(r)  =  7T(r)  (A-3) 

where  T  is  the  scalar  anomalous  potential.  For  navigation 
studies,  6g  and  T  are  the  most  likely  candidate  elements 
in  u  of  ( 3  ) .  The  disturbance  vector  represents  the 
driving  terms  in  the  inertial  navigation  velocity  error 
derivative.  The  anomalous  potential  is  related  to  geoidal 
undulation  N  by  (A-l) .  The  geoidal  undulation  could  be 
considered  an  error  term  in  a  barometric  altimeter  used  to 
stabilize  the  vertical  channel. 

For  discussion  then,  define 


(A-4) 


where  the  a-frame  is  an  arbitrary  orthogonal  frame  with 


coordinates  x,  y,  and  z. 

Using  (A-l) ,  (A-3)  can  be  rewritten  as 


\  \  \ 


Now  the  Q-matrix  function  of  (  8)  can  be  written 


0(2^ , r2)  =  $0&- T(-2^} 


where 


B2 

a2 

a2 

l 

B 

BX^BX2 

Bx-Lay2 

BXiB  z2 

g 

Bx-j^ 

a2 

a2 

a2 

1 

a 

ay jd z2 

^1^2 

ByiBz2 

g 

a*i 

a2 

B2 

B2 

l 

a 

BZiBy2 

BZfBZg 

g 

a  Zi^ 

l  a 

i  ax2 

l  a 
g  ay2 

l  a 
g  az2 

i/g2 

«(A~6) 


( (A-7) 


where  the  subscripts  on  the  x,  y,  and  z  partials  identify 
the  r  term  involved  (e.g.  S/Bx-^  means  a/ax  for  the  r^ 
terms) . 

Physical  principals  assure  that,  the  partial  deriva¬ 
tives  in  (A-5)  are  uniformly  bounded.  The  partial  deriva¬ 
tive  operations  of  (A-5)  are,  therefore,  bounded  and  linear, 
so  the  order  of  these  operations  may  be  interchanged  with 
the  linear  expectation  operator.  Then,  (A-5)  becomes 

=  3$  C0TT<£i’£2>  U  (A-8) 

1 1  2 
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where 


T(r2)]  (A-9) 

With  (A-8)  and  (A-9),  it  is  clear  that  the  Q-matrix  func¬ 
tion  can  be  produced  from  an  anomalous  potential  correla¬ 
tion  function  basis.  The  gravity  disturbance  quanties  of 
Q(r,  r ' )  are  related  to  anomaly  through  the  Stoke' s  integral 
[Ref  3:89]  and  the  Vening-Meinesz  integrals  [Ref  3;  114]]. 

So,  a  similar  line  of  reasoning  can  be  followed  to  yield 
Q(r,r')  as  a  function  of  cp  (r,r').  The  linear  operators 

00 

of  (A-8)  are  integrals  rather  than  derivatives  in  this 
case. 

With  this  demonstration  complete,  attention  can  be 
given  to  the  two  cases  which  are  likely  to  require  model 
development.  These  ares 

a.  Mission  space  geographical  restrictions  to  a 
region  known  to  have  a  residual  field  statistically 
different  from  the  global  averages,  or 

b.  Gravity  model  improvements  which  change  the  magni¬ 
tude  and  the  spectrum  of  the  residual  field  statistics. 

Even  in  these  events,  the  functional  forms  presented  in 
Section  VI  may,  yet,  be  adequate.  The  model  parameters 
would  need  to  be  re-identified  to  reflect  the  new  residual 
field  of  interest.  The  references  given  for  the  three 
functional  forms  [Refs  19,  21,  22,  and  23]]  present  some 
rationale  for  the  parameter  identification  for  the  original 
models,  and  these  suggestions  should  be  useful  in  identi- 
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fying  parameters  for  the  new  situation,  either  geographic 
or  gravity  model. 

For  the  geographically  restricted  mission,  a  geograph- 

i* 

ically  restricted,  empirical  correlation  function  could  be 
used  to  identify  the  model  parameters  on  the  basis  corre¬ 
lation  function.  Where  model  improvements  are  considered, 
the  effects  of  the  improved  gravity  model  can  be  used  to 
form  a  new  empirical  data  base.  Consequently,  a  new 
empirical  correlation  function  can  be  generated  for  subse¬ 
quent  model  parameter  identification.  Either  method 
requires  that  a  gravity  data  base  be  at  the  disposal  of 
the  analyst. 

The  required  data  may  not  be  available.  Even  if  it 
is,  the  calculations  over  all  empirical  data  would  be 
expensive.  An  alternative  might  be  considered  when  studying 
the  effects  of  a  gravity  model  improvement.  The  model 
improvement  can  be  viewed  as  a  transform  • )  over  the 
ellipsoidal  model's  residual  fields 

T '  (r)  =  J^[T(r)J  (A-10) 

where  T'  is  the  residual  field  of  the  improved  gravity 
model.  Note  thatej^(  • )  transforms  one  function  into  another 
function,  not  just  one  value  of  the  scalar  T(*)  into  another 
scalar  value.  A  suggestion  for  future  research  based  on 
this  operator  concept  is  given  in  Section  XI. 
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This  appendix  has  presented  some  of  the  statistical 
model  selection  issues.  A  categorical  answer  to  the  gravity 
disturbance  correlation  model  is  not  available,  but  several 
examples  are  developed  [Refs  19,  21,  22,  and  23]  and  are 
amenable  to  adaptation. 


APPENDIX  B 


Conversion  of  a  Spatial,  Linear  State  Space,  Statistical 
Gravity  Disturbance  Model  to  a  Temporal  Model _ 

The  time  history  of  mission  position  converts  the 
gravity  disturbance  from  a  function  of  position  into  an 
implicit  function  of  time.  The  statistical  model  for 
gravity  disturbances  undergoes  a  similar  transformation. 

For  the  linear  state  space  statistical  model,  this  trans¬ 
formation  can  be  continued  to  form  an  equivalent  temporal 
statistical  model  for  a  restricted  class  of  trajectories. 
This  conversion  to  the  time  domain  is  mentioned  in  Sections 
I  through  IV,  and  this  appendix  provides  a  derivation  to 
explain  the  assertions. 

The  linear  state  space  gravity  disturbance  model 
attempts  both  to  approximate  an  empirically  derived  corre¬ 
lation  function  and  to  yield  auto-  and  cross-correlations  of 
all  disturbance  quantities  produced  which  are  consistent 
with  gravitational  field  theory.  If,  for  example,  0*  (ljj) 

DO 

is  the  empirical  function,  we  have 

0gg(iji)  s  ^>[Ag(r1)  Ag(r2)]  (B-l) 

(rltr2)  GT 

where 

T=  C(r1,r2)Jri-r2  =  r2  cos  t)j}  (B-2) 

The  Gauss-Markov  model  is  based  on  a  state  vector  x 

-g 

which  satisfies  the  shift  invariant  relationship 

df  Cxg(t|/)3  =  Fg  xg(40  +  Gg  q( ijf)  (B-3) 
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where  q  is  a  vector  of  white  gaussian  noise  sources  satis¬ 
fying 

$Cq(i|i)  qT(4f+Ai}/)  ]  =  Qg  8(Ai|f)  (B-4) 

The  anomaly  is  an  output  of  this  process 

Ag(t|f)  =  Cgxg(<j/)  (B-5) 

The  model's  anomaly  correlation,  0” (ijr) ,  is  given  by 

OO 

0gg(A!jf)  =  CgPxg(Ai}/)  C*  (B-6) 

where 

PxgW=  <gtxg(i|r)  Xg(^A^) D  (B-7) 

which  can  be  produced  from  (B-3)  and  (B-4) .  The  difficulty 
of  performing  this  task  should  not  be  underestimated. 

The  modeler  must  form  F  ,  G  ,  Q  ,  and  C  in  order  to  pro- 

6  6  6  o 

duce  0"  (i|>).  These  model  elements  are  used  in  producing 
a  model  for  which 

0gg(iji)  ~  0gg(^)  (B-8) 

Other  disturbance  quantities  produced  by  extension  of  (B-5) 
should  have  auto-  and  cross-correlations  consistent  with 


field  theory.  A  chronicle  of  this  model’s  development  is 
provided  in  Reference  11  for  the  interested  reader. 


The  model  can  be  viewed  as 


The  error  state  transition  matrix  for  (B-l)  is 
(ty-ij;  )F 

$(i}t-i]J0)  =  e  0  g 


(B-9) 


Then, 


186 


xg(4f)  =  ®(i|>-i|f0)  2Sg(^0)  +  .fjj  $(i{J-v)Ggq(v)dv  (B-10) 
Using  the  statistical  independence  of  q(i{j)  and  x  (tjr)  with 

O 

(B-4) ,  yields  the  result 


=  »<*-*0>pg<*o>  *  <*-*0> 

+  ${l|f-v)GgQgGg  $T(t|f-v)  dv 


(B-ll) 


The  statistical  process  is  typically  modeled  to  be  at  steady 


state,  so 


pj*)  =  p_(4.0>  ■ 


(B-12) 


Since  the  derivative  of  P  with  respect  to  tjj  is  zero,  it 

O 

follows  that 


G  Q  G  =  -F  P 
g  g  g  g  g 


P  F^ 
g  g 


(B-13) 


Solution  for  Q  in  (B-13)  is  not  unique  in  general  since 
6 

G"1  could  only  exist  for  cases  where  x  and  q  have  the  same 

O  O 

number  of  elements.  The  two  Schuler  loop  cases  of  Section 
III  provide  convenient  examples  for  (B-13)  application. 

To  convert  this  spatial  model  to  the  time  domain, 
assume  i(j  changes  according  to  the  rule 


l(f(t)  =  (J;Q  +  [V(  t)/r(  t)  ]dt  =  i(j0  +  co(t)dt 


(B-14) 


where 


V(t)  is  horizontal  velocity  magnitude  and  non-negative, 
r(t)  is  position  radius  magnitude,  and 
•0)(t)  is  central  angular  velocity. 

For  this  transformation  to  be  valid  simultaneously  between 
all  points  on  the  trajectory,  the  path  has  to  lie  in  one 


.  ■  ■wi&ssit 


plane  with  a  total  change  of  less  than  180°.  Using  (B-14) 

=  [V(t)/r(t)]dt  =  o)(t)  dt  (B-15) 

Employing  the  chain  rule,  (B— 3)  can  be  converted  to 

xg[i}j(t)]  =  [u)(t)Fg]  Xgll'K*)]  +  Gg  ]£(*)  (B-16) 

where 

w(t)  =  d)(  t)q[i};{  t)]  (B-17) 

The  new  noise  correlation  is 

$  [w(t)wT(t+T)  ]  =  (i)2(  t)  Qg  fi[i];(t+T)  -  4f(t)]  (B-18) 

With  the  monotonic  assumptions  above,  the  argument  of  the 
dirac  delta  function  can  be  replaced  by  the  Taylor  series, 
first-order  approximation 
t|f(t+r)  -  ijj(t)  TU)(t) 

It  can  be  shown  that 

co(t)  5[TO)(t)3  =  8(t)  (B-19) 

So,  (B-18)  can  be  rewritten  as 

<f[w(t)  wT(t+T)]  =  a)(t)  Qg  6(t)  (B-20) 

which  completes  the  temporal  model  of  (B-l6) 


Note  that  for  V(t)=0,  x  (t)=0,  since  the  noise  strength 

O 

and  the  feedback  term  have  o)(t)  factors.  This  result  is 
intuitively  satisfying.  That  is,  when  position  is  constant, 
the  gravity  disturbance  is  constant  as  expected. 

In  order  for  the  temporal  model  to  emulate  the  corre¬ 
lations  of  the  spatial  model,  some  restrictions  must  be 
placed  on  the  type  of  trajectory  used.  Specifically,  the 


188 


trajectory  must  be  a  great  circle  of  less  than  180°,  and 
the  central  angle  must  be  monotonic  nondecreasing.  These 
trajectory  restrictions  are  further  discussed  at  the  end  of 
Section  II. 


APPENDIX  C 


Undamped  Schuler  Case,  Graphical  Results 

This  appendix  presents  the  graphical  results  of  the 
undamped  Schuler  loop  study  of  Section  III.  The  table 
below  gives  a  convenient  cross  reference  for  the  plots. 


Table  C-I 

Undamped  Schuler  Loop  Graphical  Results 
Cross  Reference 


Type  cf 

Error  Plot  Rectangular 


Figure 

Page 

Position 

Variance 

-1  *• 

192 

Pos.-Vel. 

Covariance 

C-4 

195 

Velocity 

Variance 

C-7 

199 

#  Position 
Variance 

C-10 

201 

Composite 

Position 

Variance 

c-13 

204 

Composite 

Velocity 

Variance 

C-14 

205 

Nested 

Trapezoidal _ Integrals 


Fi gure 

Page 

Figure 

Page 

C-2 

193 

C-3 

194 

C-5 

196 

C-6 

197 

C-8 

199 

C-9 

200 

C-ll 

202 

C-12 
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C-13 
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C-13 

204 

C-l*i 
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C-14 
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Figure  C-l.  Error  in  Position  Variance,  Undamped  Schuler 
Case,  Rectangular  Algorithm 
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Figure  0-4.  Error  in  Velocity  Variance,  Undamped 
Schuler  Case,  Rectangular  Algorithm 
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Figure  C-6.  Error  in  Velocity  Variance,  Undamped 
Schuler  Case,  Nested  Integrals  Algorithm 
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Figure  C-10.  Percent  Error  in  Position  Variance 
Undamped  Schuler  Case,  Rectangular  Algorithm 


Figure  C-ll.  Percent  Error  in  Position-Variance, 
Undamped  Schuler  Case,  Trapezoidal  Algorithm 


Figure  C-13-  Comparison  of  Position  Variance 
Errors  for  the  Undamped  Schuler  Case 


APPENDIX  D 


Damped  Schuler  Case,  Graphical  Results 
This  appendix  presents  the  graphical  results  for 
the  Nested  Integral  algorithm  on  the  damped  Schuler  loop 
case  discussed  in  Section  III. 
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Figure  D-l.  Error  in  Position  Variance,  Damped 

Schuler  Case 
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Figure  D-4.  Percent  Error  in  Velocity  Variance, 
Damped  Schuler  Case 
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APPENDIX  E 


Modified  Widnall-Grundy  Inertial  Navigation  Error 
Propagation  Model 

The  Widnall-Grundy  inertial  navigation  error  propaga¬ 
tion  model  [Ref  9  :26-27]  is  used  in  several  examples  in 
this  study.  A  modification  of  the  vertical  channel  feed¬ 
back  is  made  to  simulate  a  damped  vertical  channel.  This 
appendix  presents  the  specific  equations  used  to  model  the 
F-matrix. 

The  nine  error  state  vector  elements  are  defined  in 

Section  IV.  The  following  additional  terms  are  defined  for 

use  in  this  appendix. 

0  ~  design  mission  latitude 
r  ~  design  mission  radius  magnitude 
g  -  magnitude  of  gravity  vector 

e,n,z  ~  subscripts  indicate  east,  north,  or  vertical 

component  respectively 

ve>  v  ,  vz  -  earth  relative  velocity 

f  ,  ~  specific  force 

6  n  z 

(1).  earth  rotation  rate 
le 

The  following  computations  are  required  to  form  the 
F-matrix: 

^n  =  ^ie  cos  ^  (E-l) 

Qz  =  ti)ie  sin  0  (E-2) 

Components  of  local  level,  fixed  azimuth  platform 

angular  velocity  with  respect  to  the  east-north-up  frame 

are 


paq»  auwuMoy  nuqp 
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Pz  =  ve  tan  0/r 


(E-5) 


And  components  of  angular  rate  with  respect  to  inertial 


space  are 

“e  '  pe 
“n  = 

“z  =  i  *  nz 

For  later  reference,  also,  define 

kz  =  Vr 


(E-6) 

(E-7) 

(E-8) 


(E-9) 


Let  F.  .  represent  an  element  of  F.  All  elements  of  F  are 

«J 

zero  except  the  following  40  elements: 


F12  =  Pz/cos  0 
F1^  =  -  Pn/r  cos  0 

P,),  =  1/r  cos  0 


F25  =  l/r 

F36  =  1' 


P42  =  2(Qnvn  +  Qzvz)  +  Pnvn/cos  0 
F43  =  t  Pe  +  Pnkz 
*44  =  "  £  tan  0  ~  kz 

F45  =  wz  + 

F4 6  =  "wn  "  ^n 
F48  ~  "fz 
F49  =  fn 

F52  =  -ve(2nn  +  Vcos2  0> 


(E-10) 

(E-ll) 

(E-12) 

(E-13) 

(E-14) 

(E-15) 

(E-16) 

(E-17) 

(E-18) 

(E-19) 

(E-20) 

(E-21) 

(E-22) 

(E-23) 


The  "term  above  gives  the  coupling  of  6h  into  dv^.  In 
an  uncompensated  system,  this  element  is  predominated  by  a 
+2g/r  term  which  accounts  for  the  vertical  loop  instability. 
For  this  study,  altimeter  feedback  was  simply  modeled  CRef 
8]  as  a  -g/r  term.  This  modification  represents  a  stable 
vertical  loop. 

In  summary  this  appendix  provides  the  equations  that 
form  an  F-matrix  for  the  example  navigation  studies.  Equa¬ 
tions  (E-l)  through  (E-49)  supply  all  the  non-zero  elements 
for  the  Widnall-Grundy  navigation  error  state  propagation 
model  with  a  modification  to  simulate  a  damped  vertical 
channel.  The  required  inputs  to  these  equations  are  a).  , 
g,  r,  0,  ve,  vn,  vz,  fg,  fn>  and  f z.  These  data  come  pri¬ 
marily  from  the  trajectory  model  and  yield  F(t)  for  equation 
(3). 
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APPENDIX  F 


Linear  State  Space  Gravity  Disturbance  Model 


A  linear  state  space  gravity  disturbance  model  and  its 
associated  correlation  function  are  required  for  several  of 
the  examples.  An  eight-state  shaping  filter  driven  by  three 
independent  white,  gaussian  noise  sources*  has  been  developed 
CRef  193  and  gives  all  three  gravity  disturbance  terms  needed 
for  these  studies.  The  statistical  interrelationships  of 
this  model  are  consistent  with  gravitational  field  theory 
except  for  minor  approximations.  This  appendix  presents  the 
equations  required  for  both  linear  state  space  and  Nested 
Integral  covariance  analysis.  An  overall  input-output  list 
is  provided  in  the  summary. 

Figure  F-l  shows  the  shaping  filter  block  diagram  and 
output  matrix  with 


{qx(  t) 

q?(t)  . 

,3(t)J 


(F-l) 


Recall 


<£Lq(t)  qT(p)3  =  Qg( t)  8  (p-t) 


(22) 


The  noise  strengths  and  correlation  parameter  for  this  model 


were  selected  to  match  an  empirical  anomaly  correlation  func¬ 


tion  [Ref  1 93  * 


6g(t) 


g 


100 
0  4  0 

0  0  20 


(F-2) 


*  This  model  is  sometimes  referred  to  as  "third  order”  due 
to  the  three  levels  of  integration  separating  noise  sources 
from  outputs.  This  structure  forms  a  third  order  Markov  pro¬ 
cess. 
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2 

where  or  is  the  anomaly  variance  at  the  flight  altitude. 

O 

For  the  studies  which  use  this  model,  flight  altitude  and 

2 

velocity  were  constant;  therefore,  neither  a  nor  3  are 

O 

varied.  Recall, 


o  _  Y 

P  ~  A 


(85) 


where  V  is  the  horizontal  earth-relative  velocity  and  d  is 
the  correlation  parameter.  Note:  d  is  not  the  distance  at 
which  anomaly  correlation  is  l/€  of  anomaly  variance.  The 
term  correlation  parameter  will  be  used  with  the  understanding 
that  this  distance  applies  to  the  underlying  model’s  indi¬ 
vidual  feedback  gains,  not  the  resultant  filter  correlation 
distance . 


The  filter  states  in  Figure  F-l  are  named  xiQ  through 
x-j^,  since,  in  this  study,  x^  through  x^  are  navigation  error 
states  defined  in  Section  IV. 

x10(t) 

=  :  /  (f-3) 

(xi7(t)] 

The  filter  outputs  are 

T  ~  in  plane  deflection  angle, 
p  ~  transverse  deflection  angle, 

N  ~  geoidal  undulation  (not  used  here),  and 

Ag  ~  gravity  anomaly 

For  a  flight  path  heading  angle  a,  the  relationship  of 
prime,  n  ,  and  meridional,  f  ,  deflections  to  the  inplane- 
transverse  pair  is  shown  in  Figure  F-l.  Since 


+  01 


U(t)  =  jSgJ  (94) 

!6ez) 

Horizontal  elements  are  given  by  using  g  to  scale  the  deflec¬ 
tion  angles,  and  the  vertical  component  is  modeled  as  gravity 
anomaly . 


&ge  =  -g1!  =  gp  cos  a  -  gT  sin  a 
=  g  cos  a  xi;L  -  g  sin 

5^  =  -g£  =  -gp  sin  a  -  gT  cos  a 

=  -g  sin  a  x-q  -  g  cos  aUx^-gx^) 

&gz  =  -Ag  =  -g/2(  xi6-3xii4.-3xi'j7) 


(F-4) 

(F-5) 

(F-6) 


the  output  matrix 

for  (21) 

and 

(25) 

is 

g  cosa 

0  -gsina 

3g  sina 

0 

0 

0 

-g  sina 

0  -gcosa 

3g  cosa 

0 

0 

0 

(F-7) 

0 

0  0 

3g/2 

0 

-g/2 

3g/2 

gravity 

disturbance 

process 

does 

not 

begin  at  tQf 

the  filter  is  modeled  at  steady  state.  P  is  zero  at  steady 

O 

state  which  yields 


V‘>  *  - 


P 

0 

0 

0 

P 

T 

0 

0 

0 

5pr 

(F-8) 


where 


2g2  3 


(F-9) 


a2  802  43  2 

PT  =  -S  ^3  ^  3/3 

2g  L2  3/3  3/32 


(F-10) 


The  shaping  filter  state  propagation  matrices,  F  and  G  of 


(4  ),  can  be  defined  from  Figure  (F-l)  as 


°g  =‘ 


and 


where 


0 

0 

0 


0 

1 

0 


Fg  ~ 


4 


T 

0 


0  0 
0  0 
0  0 


0 

0 


T  J 


F  =  ' 

4  L  1  -3  . 


0 

0 

1 


0 

0 

0 


0 

0 

0 


g 


g 


(F-ll) 


(F-12) 


(F-13) 


and 


F  = 
T 


■-3  0  0- 

1-3  0 
.  0  1-3- 


(F-14) 


With  this  last  set  of  equations,  the  gravity  disturbance 
model  for  linear  state  space  analyses  is  complete.  The 
resulting  correlations  need  to  be  expressed  for  Nested 
Integrals.  These  correlations  will  first  be  expressed  in 
the  transverse-inplane-down  coordinate  r-frame  of  Ref  18. 
Then,  the  correlation  matrix  function  will  be  converted  by 
a  similarity  transform  to  the  east-north-vertical  n-frame 
in  which  the  navigation  analysis  is  conducted. 
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Since  the  anomaly  correlation  is  isotropic  in  heading, 
cross  correlations  of  transverse  deflection  with  the  other 
terms  are  zero.  This  result  can  be  deduced  from  Figure  F-l 
where  (j.  is  produced  by  a  statistically  independent  shaping 
filter.  The  Q-function  has  the  form  for  r-^  and  rv,  sepa¬ 
rated  by  central  angle  I 


QT(£!»r2)  = 


0 


0 


(ijr) 

0  0W  0W 

“  “  rrrv- 


TT 


er 


0  0W  0W 

Tg  ^gg  J 


(F-15) 


where  the  superscript  r  indicates  the  coordinate  frame  and 
where,  for  example, 


4S1  ■  <£Cn(r)  h<s’>] 


(F-16) 


the  expectation  is  over  all  (r,r’>  pairs  in  the  sample  space 
which  are  separated  by  i{j  central  angle. 

For  brevity  define 


M  = 

Then,  2 

=  Zfi-  (i  +  m)  e"M 

^  2g2 

2 

0(^)  =  -&X  (1  +  M  -  M2)  e"M 

TT  2g2 

0gg}  =  ag  (1  +  M  “ =m2)  e'M 

+  M) 

=  -  0W 

gT  % 


(116) 

(F-17) 

( F-18) 

(115) 

(P-19) 

(F-20) 


2_2 


0(tJ/)  =  0(!j/)  =  =  -  o  (F-21) 

vT[l  M-T  v&  \i&  '  ' 

Using  (116)  through  (F-21),  Q(r,rj)  can  be  formed  with  (F-15)  . 
This  correlation  function  satisfies  a  need  of  Nested  Integral 
analysis.  Qtr-^rg)  is  formed  in  the  T-frame,  inplane- 

transverse-down,  coordination  and  the  navigation  error  terms 

» 

are  in  an  n-frame,  east-north-vertical.  An  arbitrary  vec¬ 
tor  Z,  expressed  in  n-frame  coordinates  is  written  Zn. 

Define  C?  as  the  coordinate  transformation  matrix.  Then 
m 

Z11  =  ZT  (F-22) 


Recall  that  heading  angle  a  separates  the  north  and  inplane 
axes  in  the  horizontal  plane,  but  note  that  a,  here,  is  not 


is  given  by,  paraphrasing  (F-22), 


east 

-cos  a  sin  a  0 

transverse 

north 

= 

sin  a  cos  a  0 

inplane 

vertical 

I - 

O 

O 

1 

M 

1 _ 

down 

Since, 

C. 


=  c  =  c  ^ 

nr  m  Tn 


(F-23) 


(F-24) 
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(P-25) 


Qn(£l-£2)  -  CnT  «T(£l,r2)  CnrT 

This  Q(r^,r2) ,  then,  is  the  form  required  to  interface  cor¬ 
rectly  with  the  Widnall-Grundy  navigation  error  model.  The 
coordinate  transformation  performed  here  could  have  been 
performed  in  the  G-matrix  with  a  redefinition  of  u  in  (110) 
to  the  t- frame  rather  than  the  present  n-frame  coordinates. 

The  heading  angle,  a,  of  (F-23)  must  correspond  to  the 
navigation  coordinates  in  effect  at  that  time.  The  argu¬ 
ments  of  Q(rltr2)  are  two  positions,  rtt-^)  and  r(t2).  The 
great  circle  arc  from  r(t-^)  to  r(t2)  has  heading  at  the 
r(  t-^)  end  and  a2  at  the  r(t2)  end.  In  the  Nested  Integrals 
approach,  for  t^  <  tn,  Q(tn,  t^)  evaluations  are  r(tn)  and 
rCt^).  an  must  be  used  since  t=tn  when  these  calculations 
are  performed  in  these  evaluations.  If  position  coordinates 
are  given  in  longitude  latitude  0,  and  altitude  h,  the 
heading  angle  needed  is 

_1  j  sin  Un-Xi)  cos  0i 


*n 


=  tan 


cos  0^  sin  0n  cos  ( Xn~X^)  -  sin  0^  cos  0n 

(P-26) 


The  central  angle  between  r(t^)  and  r(tfi)  is 

tj/^n  =  cos-1  Ceos  0^  cos  0n  cos  (Xn-X^)  -  sin  0^  sin  0  ] 

(F-2?) 

These  complete  the  requirements  for  Nested  Integrals  analysi 
In  summary,  this  appendix  has  presented  a  previously 
developed  [Kef  1 9]  eight-state,  linear  shaping  filter  model 
for  the  gravity  disturbance  process.  For  the  linear  state 
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o 

space  covariance  analysis  method,  inputs  are  Vf  d,  g,  a*9 

D 

and  a.  The  outputs  are 

Fg(t)  of  (  4  )  and  (9  )  given  by  (F-12)  through  (F-14) 
and  (85) 

6  (t)  of  ( 4  )  and  (10)  given  by  (F-ll) 

O 

C  (t)  of  ( 5  )  and  ( 9  )  given  by  (F-?) 

0 

Q  (t)  of  (  6  )  and  (14)  given  by  (F2) 

Part  of  the  initial  condition  for  (14)  is  given  as 
P  (tQ)  in  (F-8)  through  (F-10). 

p 

The  Nested  Integrals  method  inputs  are  g,  d,  o  ,  cos  0  , 

D  ** 

cos  0^,  sin  0n,  sin  0^,  and  Xn,*  the  output  is 

Q(rn,r^)  in  n-frame  coordinates  using  (115) ,  (116), 
and  (F-l?)  through  (F-25)  where 

central  angle  tjj^n  is  given  by  (F-27)  and 

arc  heading  an  is  given  by  (F-26) . 


APPENDIX  G 


Great  Circle  Flight  Path 

The  trajectory  generator  is  required  to  yield  position, 
velocity  and  specific  force  throughout  the  interval  of  study. 
A  general  trajectory  generation  program  [e.g.,  Ref  ?  ]  could 
be  used  for  this  purposes  however,  computer  resources  can 
be  conserved  for  this  study  by  programming  the  following 
closed-form  solutions.  An  input-output  list  is  provided  in 
the  summary  to  this  appendix. 

A  spherical  earth  of  radius  Rg  is  assumed.  Two  sets  of 
earth-centered,  earth-fixed  coordinates  facilitate  the  devel¬ 
opment.  The  a-frame  is  defined  by  x-axis  through  the  equator 
where  the  great  circle  crosses  while  the  path  is  moving 
north.  The  z-axis  is  along  the  north  polar  axis  (parallel  to 
£iiie)  ,  and  the  y-axis,  in  the  equatorial  plane,  completes  the 
right-handed  coordinate  set.  The  b-frame  has  the  same  x- 
axis  but  the  y-axis  lies  in  the  great  circle  plane.  Figure 
G-l  shows  this  basic  geometry. 

The  basic  angles  of  interest  in  Figure  G-l  are 
i  -  the  earth  relative  inclination  of  the  path 
\  -  the  longitude  from  x0 

cx 

<t>  -  geocentric  latitude 
A-  great  circle,  earth  central  angle 
For  brevity,  time  arguments  are  not  given.  The  radius  magni¬ 
tude,  r,  at  constant  altitude,  h,  is 

r  =  Re  +  h  (G-l) 


HaCEDIWQ  PAQB  BUMmHOT  HUfl) 
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Figure  G-l.  Great  Circle  Geometry:  Earth-Centered, 
Earth-Fixed  Coordinates 


Since  h  is  constant,  r  is  constant.  Now  since  V  and  h  are 
constant 

A  =  F  +  Ao  (G"2) 

The  central  angle  tji(t2>t^)  between  r(t^)  and  r(t2)  is  simply 

liKtg.tj)  =  |A(t2)  -  A(t1)|  =  -  I  t2  -  txl  (G-3) 
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:!! 

t 

H  i 


(G-4) 


The  coordinate  transformation  from  b-frame  to  a-frame  is 

0  0 


'ba 


1 

0 

0 


cos  1 
sin  i 


-sin  l 
cos  i 


(G-5) 


So, 


lCosA 

r“  =  r  <*cos  i  sinA, 
sin  i  sinA 

-1 


X  =  tan-i  (yQ/x_)  =  tan-'*"  (cos  i  tanA) 

a  a 

cos  0  =Vxf  +  =  V 


2  2.  .  2 
cos  +  cos  l  sin 


sin  0  =  —  =  sin  i  sin  A 


(G-6) 

(G-7) 

(G-8) 

(G-9) 


From  Figure  G-2,  the  heading  angle  a  is  defined.  For  a 
great  circle 


cos  i  =  sin  a  cos  0 
sin  a  =  cos  i/cos  0 


(G-10) 

(G-ll) 


-sm2a 


cos  a  =< 


(G-12) 


-90°  <  A  <  90° 

1-^  l-sin2a  ^  90°  <  A  <  2?0° 

Since  the  velocity  vector  is  of  magnitude  V  and  lies  a 
radians  clockwise  from  the  north  axis  in  the  horizontal 
plane , 

vn  =  V  cos  a  (G-13) 

and 

ve  =  V  sin  a  (G-14) 


1  V  =  earth  relative  velocity 

i  -  subscript  means  with  respect  to  inertial  frame 

e  -  subscript  means  with  respect  to  earth-fixed  frame 

f  -  is  specific  force  vector 

g  -  is  gravity  vector 

(1).  _  -  is  earth  rotation  vector 
“ie  - 


Using 


I  =  Y  +  2(oie  XV  -  g 


(G-17) 

(G-18) 

(G-19) 


-2Vu)ie 

sin  0 

cos 

a  1 

2V^ie 

sin  0 

sin 

a  < 

l  (G-20) 

V* 

e-~ 

-2VWie 

cos 

‘J 

. 

Equation  (G-l) ,  (G-7),  (G-8)  ,  (G-9)  provide  the  position 
coordinates  needed  for  Nested  Integrals  correlation  evalu¬ 
ation.  Latitude  0  need  not  be  found  directly  since  only 
sin  0  and  cos  0  are  required  for  the  i{;  and  calculations. 
For  F-matrix  evaluation,  these  equations  are  supplemented 
by  the  relative  velocity,  (G-13)  through  (G-15) ,  and  spe¬ 
cific  force  (G-20) . 

In  summary,  a  great  circle  trajectory  has  been  modeled 
in  sufficient  detail  to  support  the  covariance  analyses  of 
this  study.  Inputs  are  t,  t0  R  ,  h,  V,  A_,  i,  and  g. 

”  c  U  -L  c 
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Outputs  are  r,  X,  cos  0,  sin  0,  Vn,  sin  a»  cos  a,  and  f”. 
These  outputs  interface  with  other  models  and  analyses  as 
follows: 

1)  g»  r,  cos  0,  sin  0,  Vn,  and  f11  for  F-matrix  evalu¬ 
ation  of  Appendix  E 

2)  cos  a  and  sin  a  for  C(t)  evaluation  of  (F-7)  in 
Appendix  F 

3)  X,  cos  0  and  sin  0  for  Nested  Integrals  position 
data  to  be  used  in  and  ij/^n  calculations  of 
(F-26)  and  (F-2?)  in  Appendix  F. 

4)  h  for  r2)  evaluations  for  the  Tscherning-Rapp, 

and  the  Heller-Thomas  correlation  models  used  in 
Sections  V  and  VI  (see  also  Appendices  J  and  K) . 
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APPENDIX  H 

MINOR  CIRCLE  FLIGHT  PATH 

The  minor  circle  flight  path  trajectory  model  is 
presented  in  this  appendix.  Based  on  input  data  which  shall 
be  summarized  later,  this  model  provides  the  position, 
velocity  and  specific  force  information  required  for  the 
covariance  analyses  used  in  this  study. 

The  basic  geometry  and  the  earth-centered,  earth-fixed 
coordinates  (a-frame)  are  pictured  in  Figure  H-l.  In  this 
figure 

\  is  longitude  from  x_ 

CL 

0  is  geocentric  latitude 

A  is  minor  circle  angle 

r  is  radius  vector 

r  is  minor  circle  radial  vector 
— c 

V  earth  relative  velocity 

r^  fixed  radial  vector  to  center  of  minor  circle 
The  radius  magnitude  is  constant 

r  =  Re  +  h  (H-l) 

where  Rg  is  earth  radius  and  h  is  the  constant  altitude. 

The  minor  circle  constant  angular  velocity  is 

W  =  J-  (H-2) 

c 

Now  the  minor  circle  angle  is  given  by 

A  =  0)(t-to)  +  Aq  (H-3) 

where  t  is  initial  time  and  AQ  is  initial  angle.  Note  that 
time  arguments  are  dropped  for  brevity. 
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Now  given 


(H-4) 


(H-5) 


Also, 


r  °  i 

=  r  J  sin  A  ► 
— c  c  ] 

l^cosA^ 


(H-6) 


These  results  complete  the  minor  circle  trajectory.  The 
required  inputs  are  t,  tQ,  Re,  h,  rc,  V,  AQ,  g,  and  wie. 

The  outputs  are  as  follows*  g,  r,  cos  0,  sin  0,  vn,  and 
fn  for  F-matrix  calculations;  a  for  C(t)  calculation  in 
(G-7)  for  linear  state  space  covariance  analysis;  and  cos  0, 
sin  0,  h,  and  X  for  Nested  Integrals. 
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The  equation, (113),  for  central  angle  presented  in 


Section  IV  can  be  derived  by  using 


^1,2  =  cos’1^£1  •  £2^rlr2^  (H-21) 

where  r-^  is  rCt-^)  and  r2  is  r(t2) 

when  (H-7)  is  used  for  both  r^  and  r2,  the  result  is 


^1,2  " 


=  cos 


-1, 


r  - 


r^Cl-cosCp-  ( t2- t-j^)  }3 
_ c _ 


(H-22) 
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Circular-Error-Probable.  Computation 

The  analyses  performed  in  Sections  IV  through  IX  used 
a  nine-state  navigation  error  model  and,  therefore,  produced 
an  8l-element  covariance  matrix.  Since  this  matrix  is 
symmetric,  only  45  elements  need  to  be  considered  in  any 
comparative  study.  In  the  verification  study  of  Section  IV, 
all  45  independent  elements  are  presented  for  the  great 
circle  case  comparison.  Such  an  exhaustive  comparison  is 
rarely  necessary  since  the  primary  interest  is  in  position, 
and  perhaps  velocity,  statistics.  The  desire  for  a  more 
compact  performance  index  can  be  met  with  some  second  order 
statistic  of  position. 

Several  alternatives  are  available,  but  circular-error- 
probable  (CEP)  is  the  most  common.  This  statistic  gives  the 
horizontal-plane  radius  which,  when  centered  on  the  mean, 
encompasses  one-half  the  population.  This  appendix  presents 
the  method  used  to  calculate  CEP  in  the  studies  of  Sections 
IV  through  IX . 

CEP  can  be  calculated  by  two  meth  ds.  The  entire  popu¬ 
lation  can  be  accounted  for  by  the  instructions  above,  or 
the  frequency  function  can  be  integrated  until  the  radius 
is  found,  which  gives  a  probability  of  one-half.  Neither 
the  population  nor  the  frequency  function  are  known  in  gene¬ 
ral.  Since  the  navigation  error  propagation  is  a  linear 
process,  the  cases  which  employed  the  gaussian  noise  models 
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will  have  multidimensional  gaussian  distribution  for  the 
error  states.  In  these  instances,  a  rather  straightforward 
method  is  available  to  estimate  CEP. 

The  population  is  not  known  for  the  other  cases,  but 
the  errors  which  drive  the  navigation  covariance  are  numerous. 
That  is,  the  distributed  earth  mass  nonuniformity  can  be 
construed  as  millions  of  error  sources,  and  the  gravity 
disturbances  are  the  sum  of  these  errors.  Appealing  to 
the  Central  Limit  Theorem  [Ref  15:109],  the  error  population 
is  probably  well-approximated  by  a  gaussian  distribution. 

If  the  gaussian  distribution  is  used  as  an  approximation 
for  cases  of  unknown  population,  then  the  same  method  can  be 
used  to  calculate  CEP  in  all  instances. 

Of  the  45  independent  elements  of  the  navigation  error 
covariance  matrix,  only  three  are  required  in  the  CEP  calcu¬ 
lation.  These  elements  are  the  latitude  and  longitude  error 
covariances.  The  ultimate  CEP  answer  will  be  in  linear 


measure,  so  these  angular  errors  must  be  converted  to  linear. 
The  east  component  of  position  error  is 

6xg  =  r  cos  0  6X.  (1-1) 

and  the  north  component  is 


(1-2) 

(1-3) 


and 


24  0 


r2  cos20  <g*(  6X2) 
r2  cos  0  60) 


r2  cos  0<f(6A.  60) 
r2.  <f(602) 


(1-4) 


The  three  expectations  needed  to  evaluate  (1-4)  come  directly 
from  the  navigation  error  covariance  matrix  for  the  state 
vector  definition  used  in  Sections  IV  through  IX. 

For  the  bivariate,  normal  distribution  described  by 
Pp  from  (1-4),  CEP  was  approximated  using  the  following 
procedure  CRef  20  ]s 
Define, 


a2  =  $(6x2)  =  r2  cos20  <f(6\2)  (1-5) 

a2  =  (gUx2)  =  r2  <f(602)  (1-6) 

-  5V  &*n>  -  ^  <f(6x  <J-7) 

v  XI  v  XX 

A  new  coordinate  frame,  x-y,  is  defined  by  rotating  y  from 
the  east-north  frame  until  the  correlation  coefficient  is 
zeros 


For  this  condition, 

45° 

I 

Y  = 


tan 


-1 


ae  "  CTn 


2  ae  an 
2  2 
°e  “  °n  J 


(1-8) 


(1-9) 


ae  ^  °n 
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Then, 


= 

4 

2 

cos  y 

+  0e0n  sin 

2y  + 

2  2 

0„  sin^Y 

(1-10) 

4’ 

4 

•  2 
sin  v 

-  °ean  sin 

2Y  + 

2  2 
an  cos  y 

(1-11) 

Next,  the  maximum  and  minimum  variances  are  identified. 

amax  =  Maximum  tox,  ay}  (1-12) 

°min  =  Minimum  Cax,  oy]  (1-13) 

Finally  the  estimate  is  made  by 

CEP  =  0.562  +  0.615  0Bin  (1-14) 

In  essence,  (1-5)  through  (1-7)  transform  the  statistic 
to  linear  measure.  Next,  (1-9)  gives  the  resolvent  angle 
with  which  prime  variances  are  calculated  by  (1-10)  and 
( I— 11) .  Then,  the  prime  variances  are  ordered  by  (1-12) 
and  (1-13),  and,  finally,  the  CEP  calculation  performed  by 
(1-14).  The  input  data  for  these  calculations  are  r,  cos  0, 
$<.  6i2).  (rf(  60^) ,  and  ^(6\  60).  The  output  of  these  compu¬ 
tations  is  an  estimate  of  circular-error-probable. 

An  exception  to  the  above  rule  is  the  CEP  calculation 
for  the  Nested  Integrals  results  of  Section  V.  The  Monte 
Carlo  results  cited  in  Section  V  were  based  on  downrange  and 
crossrange  variances.  The  above  procedure  can  be  used  to 
yield  the  desired  comparable  result  if  ( I-9)  is  not  applied. 
For  a  downrange  -  crossrange  resolution  of  the  data,  simply 
use  the  complement  of  the  heading  angle  in  place  of  the  (1-9) 

result.  For  Section  V,  then,  the  heading  angle  is  an  addi¬ 
tional  required  input. 


APPENDIX  J 


Anomaly  Decree  Variance  Correlation  Model 

The  statistical  analysis  of  navigation  errors  induced 
by  gravity  disturbances  can  be  based  on  a  disturbance  corre¬ 
lation  model  as  discussed  in  Sections  I  and  II.  Correlation 
models  are  also  used  in  smoothing  and  predicting  problems, 
and  geodesists  have  used  this  approach  in  processing  gravity 
measurements.  These  gravity  measurements  come  in  many  forms 
(e.g.,  anomaly  and  vertical  deflections),  and  all  measure¬ 
ments  cannot  be  made  on  a  convenient  surface.  Therefore, 
some  advanced  statistical  models  for  geodetic  data  proces¬ 
sing  have  the  capability  to  produce  all  the  needed  correla¬ 
tions  for  navigation  analysis. 

The  Tscherning-Rapp  anomaly  degree  variance  model  CRef 
2i]  was  developed  for  a  geodetic  application.  This  corre¬ 
lation  model  was  used  in  the  Monte  Carlo  analysis  [Ref  10J 
and  is  used  in  Sections  V,  VI  and  VII  of  this  work.  This 
appendix  provides  some  background  on  this  model  and  on  the 
manner  in  which  it  was  used  in  these  navigation  analyses. 

The  anomaly  correlation  function  from  the  geodesist's 
viewpoint  is  the  mean  over  some  region  of  anomaly  products. 
The  ensemble  of  missions  which  gives  the  sample  space  basis 
for  inertial  navigation  studies  is  not  the  geodesist's  per¬ 
spective.  To  him 

cpcc(r,r')  =t/tf[tig(r)  Ag(r’)] 

(r,r'MT 


(J-l) 


where 


Jt  is  the  mean  operator  over  the  set T which  restricts 
(r,r')  pairs  (e.g.,  to  a  (j;  central  angle  separation). 

For  continuous  data  over  a  spherical  surface  and  central 


angle  as  the  shift  variable, 


T  =  {(r,r')|r*r'  =  r^  cos  t|/} 


(J-2) 


and 


*«<♦>  *  h  55  *«<?•*> 


X=0  0=-w/2 


a-0  cos  0  d0  dX  da 


(J-3) 


where 

0  is  geocentric  latitude 

X  is  longitude 

a  is  heading  at  r  of  the  great  circle  to  r’ . 

This  difference  in  approach  in  defining  the  correlation  is 
artificial.  In  fact,  the  statistics  defined  by  (J-3)  are 
precisely  the  sort  which  "represent"  the  anomalous  field 
over  which  the  navigation  mission  occurs. 

The  form  of  (J-3)  does  define  a  function  with  both 
homogeneous  and  isotropic  features  (see  Appendix  A).  The 
integrals  over  all  0  and  X  produce  a  global  average  consis¬ 
tent  with  a  homogeneous  assumption.  The  integral  over  all 
a  gives  a  function  which  represents  an  average  for  all 
headings  -  isotropic  assumption  [Ref  4s 853- 

Global  data  is  required  to  approximate  the  integrals  in 
(J-3) •  Reference  21  discusses  the  problems  encountered  in 
estimating  qp  (ijj)  from  a  restricted  data  base,  but  empirical 

oo 

correlation  functions  have  been  produced  in  spite  of  diffi- 
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culties.  The  modeling  approach  is  to  create  a  closed-form 
mathematical  model  which  is  capable  of  approximating  this 
empirical  correlation  function.  Then,  one  derives  other 
needed  auto-  and  cross-correlations  from  this  basis  (see 
Appendix  A) . 

The  correlation  function  of  (J-3)  can  be  expanded  into 
a  series  of  (zonal  spherical  harmonics)  Legendre  polynomials 
as  CRef  4:85] 

00 

cp^(ilj)  =  2  cn  Pn  (cos  ^  (J-*0 

ss  n=0  n  n 

where 


P  ( • )  is  the  Legendre  polynomial 
cn  is  the  associated  coefficient. 

The  name  "anomaly  degree  variance"  comes  from  the  observa¬ 
tion  that  each  cn  represents  a  contribution  to  the  anomaly 
variance.  This  form  can  be  generalized  to  points  off  the 
reference  surface  by  defining 


s  = 


si 

rr 


(J-5) 


where 

R  is  the  spherical  surface  used  as  a  datum  (Bjerhammer 
sphere  J^Ref  21:52]). 

Whereas  anomaly  is  not  harmonic  (satisfies  LaPlace* s 
Equation),  the  product  r*Ag  is  harmonic  CRef  3:90],  So, 
the  general  form  for  (J-4)  is 

(W<I')  =  JQ  cn  ®n+2  Pn(cos  ^  <J'6) 
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1 

1 


Note  that  altitude  information  enters  (J-6)  through  the 
radius  magnitudes  in  (J-5)*  The  statistics  of  (J-6)  are 
therefore  nonstationary  with  respect  to  altitude  change,  but 
for  a  fixed  altitude,  the  statistics  are  stationary. 

Low  frequency  coefficients  of  (J-6)  can  be  identified 
using  the  empirical  data  base.  The  high  frequency  terms 
must  at  some  point  be  approximated.  All  cn  are  non-negative 
[Ref  4s 85]  and  some  are  positive  giving  a  positive  definite 
variance 


2 

n=0 


n 


( J-7) 


This  equation  demonstrates  that  cfi  represents  the  anomaly 
variance  contribution  from  the  n-th  degree. 

One  approach  to  modeling  the  anomaly  correlation  is  to 
create  a  mathematical  form  which  will  fit  the  low  ordered 
cn's  and  prescribe  the  manner  in  which  high  frequency 
coefficients  approach  zero.  Tscherning  and  Rapp  propose 
five  such  rules  in  Reference  21  and  develop  models  for  four 
of  these.  These  complete  models  include  auto-  and  cross¬ 
correlations  for  anomaly,  deflection  angles,  and  geoidal 
undulation.  Each  of  these  models  is  based  on  an  algebraic 
expression  relating  cn  to  n. 

The  fourth  rule  was  completely  developed  and  program¬ 


med  for  computer  application.  That  model  is,  also,  recom¬ 
mended  in  the  conclusion  [Ref  21s30] 


8 


cn  = 


A(n-l) 
[n-2) (n+B 


for  n  <  2 
for  n  >  2 


(J-8) 


where  A  and  B  are  constants  identified  to  approximate  the 
empirical  correlation.  Then, 

cp„J<l>)  =  2  ( —  T+2  P.(cos  A)  (J-9) 


■  j3  '(ri-M'J  ( Jr)"*2  pn<cos  *> 


The  values  for  A,  B,  and  R  were  determined  from  empirical 
data  and  are  presented  in  Table  J-I,  below. 

Table  J-I 

Model  4  Data 

Symbol _ Value _ Units 

A  425.28  mgal2 

B  24  (exact) 

R  6369.8  *  km 

*  based  on  a  637I.O  km  mean 
earth  radius 

Equation  (J-9)  is  obviously  not  a  simple,  computation¬ 
ally  efficient,  closed-form  mathematical  expression.  One 
property  of  Legendre  polynomials  redeems  this  modeling 
choice.  Consider  the  triangle  formed  by  r  and  r'j 


J 


By  the  law  of  cosines 

a2  =  r2  +  r'2  -  2rr'  cos  ij; 


1 

a 


■V 


(J-10) 

(J-ii) 


r%  1  -  2bs  +  s 


-  2rr '  +  r’^ 
where,  for  now, 

s  =  r'/r  ( J-12) 

and 

b  =  cos  4/  (J-13) 

Expanding  the  right  hand  side  of  (J-ll)  as  a  power  series  in 
s  gives  [Ref  302-33] 


Vl  -  2bs  +  s2 


=  l  2  sn  P  (b) 
r  n=0 


=r  =  2  Sn  P  (cos  d;) 

.2  n=0  n 


V 1  -  2s  cos  ijj  +  s‘ 

Now  (J-9)  can  be  written  as 

cp  (tjj)  =  s2  A  2  — -  sn  P  (cos  i|j) 

where 


(J-l4a) 


( J-I4b) 


n=0  (n-2) (n+B) 


s  = 


R 


rr 


By  partial  fraction  expansion 

.nzl _ _  =  lABj-21  (B+l)/(B+2) 

(n-2) (n+B)  n-2  n+B 

So  a  closed-form  expression  for  cp  can  be  found  if 

J0  "wT  Pn(cos  ^ 


(J-15) 


(J-16) 


(J-17) 
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can  be  given  in  closed  form  using  ( J-l4b) .  Expressions 
for  these  terms  are  derived  in  Reference  21  by  using  the 
fact  that 

t  n+i 

So  sn  ds  =  s^iT-  ,  n+i  >  0  ( J-18) 

The  derivations  are  presented  in  Reference  21,  pages  30 
through  46,  and  will  not  be  repeated  here.  The  above 
discussion  shows  the  general  approach  taken. 

A  computer  subroutine  COVA  was  developed  and  documented 
in  Reference  21,  pages  85  through  89.  This  subroutine  uses 
Model  4  and  the  data  in  Table  J-I  to  produce  the  following 
correlations  as  a  function  of  central  angles 

a.  Anomaly  autocorrelation, 

b.  Anomaly  and  downrange  deflection  cross-correlation, 

c.  Anomaly  and  geoidal  height  cross-correlation, 

d.  Downrange  deflection  autocorrelation, 

e.  Crossrange  deflection  autocorrelation, 

f.  Downrange  deflection  and  geoidal  height  cross¬ 
correlation,  and 

g.  Geoidal  height  autocorrelation. 

Of  these  a,  b,  d,  and  e  are  used  in  Nested  Integrals 
analyses  of  Sections  V,  VI,  and  VII. 

Since  the  intent  of  this  model  is  to  provide  high 
frequency  cn's,  provision  was  made  to  specify  up  to  the 
first  300  coefficients.  The  closed  form  expression  for 
the  infinite  summation  is  altered  by  removing  the  effects 
of  each  coefficient  to  be  replaced  and  adding  in  the 


effects  of  the  user-supplied  coefficients.  This  feature 
is  used  in  the  spherical  harmonic  modeling  study  of  Section 
VII. 

In  performing  the  Monte  Carlo  study  discussed  in 
Section  V,  Geodynamics  discovered  minor  errors  in  COVA 
[Ref  10tl83  and  these  changes  should  be  made  before  the 
subroutine  is  used  on  a  trajectory  with  altitude  variations. 
The  subroutine  used  in  this  study  was,  in  fact,  programmed 
by  Geodynamics.  The  subroutine  was  verified  by  recreating 
Table  10  of  Reference  21. 

The  anomaly  degree  variance  correlation  model  developed 
by  Tscherning  and  Rapp  offers  an  alternate  form  of  statis¬ 
tical  model  for  use  in  the  analysis  of  navigation  errors 
induced  by  gravity  disturbances.  The  corrected  subroutine 
COVA  offers  a  convenient  useful  tool  for  direct  use  of  this 
model  in  gravity  model  performance  evaluations. 
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APPENDIX  K 


Attenuated  White  Noise  Correlation  Model 

For  long  range  missions  such  as  intercontinental 
ballistic  missiles,  the  effects  of  earth  curvature  must  be 
considered  in  navigation  performance  analyses.  The  linear 
state  space  gravity  disturbance  correlation  model  is  based 
on  a  tangent-plane,  flat-earth  approximation.  Such  flat- 
earth  models  are  of  questionable  validity  on  long-range 
mission  analysis.  The  attenuated  white  noise  correlation 
model  was  developed  to  provide  a  round-earth  alternative. 
This  appendix  provides  a  partial  summary  of  this  develop¬ 
ment  and  lists  the  correlation  functions  used  in  the  attenu¬ 
ated  white  noise  case  of  Section  VI. 

The  gravity  disturbance  process  on  the  earth's  surface 
can  be  statistically  summarized  in  the  form  of  the  anomalous 
potential  correlation  function  (see  Appendix  A) . 

t(e')3  (K-l) 

The  anomalous  potential  is  a  spatially  correlated  process. 

A  ploy  used  in  modeling  a  temporally  correlated  process  is 
to  view  the  random  occurrences  as  the  output  of  a  linear 
system  driven  by  white  gaussian  noises.  Symbolically, 


white,  temporal  shaping  correlated,  temporal 

noise  filter  output 
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The  approach  used  by  Heller  [Ref  23]  in  producing  the  attenu¬ 
ated  white  noise  was  to  model  the  earth-surface  anomalous 
potential  as  the  upward  continuation  of  a  white  noise  pro¬ 
cess  over  a  subterranean  sphere  CRef  4:4-5].  Symbolically, 


T(P) 


* 


Upward 

Continuation 

Integral 


T(R) 


white,  spatial 
noise 


Laplace's  correlated,  spatial 
Equation  output 


where 


P  is  radius  vector  on  the  sphere,  and 
R  is  radius  vector  on  the  earth's  surface. 


constant  variance  given  by 

%.(£.£’)  =  0TT(i|f)  =  A  f (({;)  (K-2) 

and  f(i{>)  is  the  unit  impulse  function  on  the  sphere  satis¬ 
fying 
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(K-3) 


2 irp2,  f(i|j)  sin  i{/  dt|/  =  1 

or 


f(i|»)  =  - -  (K-4) 

2 irp  sin 

The  anomalous  potential  above  the  noise  source  is  given  by 
upward  continuing  (K-2) ,  An  extensive  derivation  is  given 
in  Appendices  A  through  F  of  Reference  23-  The  result  is 
0T.r<r.£')  ~  0TT(rr'f!j/)  = 

d^R-dPUrr')2  -  (R-d)4]  g2 _ 

CR4-(R-d)4]C(R-d)4+(rr,)2-2rr,(R-d)2  cos  ^]3/2  (K~5) 

p 

where  g£  is  the  earth  surface  variance  level  (let  r=r’=R 
and  ij;=0 )  .  A  complete  set  of  disturbance  terms  and  auto- 
and  cross-correlations  were  derived  from  (K-5)  by  methods 
outlined  in  Appendix  A  of  this  work.  The  key  benefit  of 
this  model  is  that  (K-5)  correctly  accounts  for  altitude 
changes  in  r  or  r'  without  requiring  a  new  integral  evalu¬ 
ation  [Ref  23:4-6]. 

The  depth  parameter  d  acts  in  a  similar  manner  to  the 
feedback  gain  in  the  linear  state  space  model.  Correlation 
distance  is  directly  related  to  d.  For  a  d  of  10  nm,  the 
earth  surface  arc  correlation  distance  is  11.15  nm. 

When  d  is  500  nm,  correlation  distance  is  468  nm  [Ref  23: 
pF-13]. 

Asymptotic  forms  of  these  round-earth  correlation  models 
were  also  developed  in  Reference  23.  These  asymptotic  forms 
come  from  first  order  expansions  in  i|j  and  in  d/R  [Ref  23» 


■  ■  f-  (  r  i- 
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4-9  through  4-11],  These  mathematically  simpler  expressions 
were  used  in  the  Section  VI  study. 

The  asymptotic  version  of  (K-5)  is 


0M(r,r'4)  ■  .  ^,2' W') 


C(2d+h+h')2  +  R2  i|>2]3^2  T 


(K-6) 


where 


h  =  r-R 


and 

h  =  r  -  R. 

For  this  study  0TT  was  not  required.  The  associated  model 

correlations  which  were  used  are  as  follows: 

2 

Define  a  the  earth  surface  gravity  anomaly  variance  value 

o 


as 


Let, 


-3  2 

2  _  3CTT 

0CT  ~  2 

g  2d^ 


T  =  Rijj 

A  =  2d  +  h  +  h ' 

B  =  8d4o|/(A2+T2)3//2 
C  =  A  •  B 

Then  the  transverse  disturbance  autocorrelation  is 
0|i(1(r,r’4>  =  C(A2+t2) 

The  inplane  disturbance  autocorrelation  is 
0TT(r,r’  ,i|i)  =  C(A2-4t2) 


(K-7) 

(K-8) 

(K-9) 

(K-10) 

(K-ll) 

(K-12) 

(K-13) 


The  inplane  and  radial  disturbances  cross-correlation  is 
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0Tg(r(r'4)  =  tB(4A2-t2)  (K-14) 

The  radial  disturbance  autocorrelation  is 

0gg(r,r*4)  =  C(2A2-3t2)  (K-15) 

p 

Given  (K-15)#  the  modeler  has  two  parameters,  a  and 

D 

d,  to  use  in  approximating  an  empirical  correlation  func¬ 
tion.  A  more  flexible  model  can  be  created,  however,  by 
simply  viewing  the  white  noise  shell  which  gave  (K-15)  as 
one  of  a  set  of  statistically  independent  processes.  With 

a  ct  and  a  d-  for  each  of  n  shells,  the  modeler  has  2n 
gi  1 

parameters  with  which  to  build  an  overall  correlation  model: 
For  example, 

n 

0-Ar,r'  ,ijj)  =  2  <S_._  (r,r*,\|f)  (K-16) 

gg  i=1  ggA 

The  other  auto-  and  cross-correlations  are  summed  in  the 
same  manner.  A  three-level  model  was  proposed  and  para¬ 
meters  identified  in  the  Reference  23  study.  These  para¬ 
meters  were  used  in  the  Section  VI  example  and  are  given 
in  Table  K-I. 


Table 

k-1 

Parameters  for  the  Three-Level 
Attenuated  White  Noise  Gravity 
Disturbance  Correlation  Model 

i 

di 

n.m. 

% 

mgal2 

Percent  of 

Total  Variance 

1 

1002 

175 

9.6 

2 

179 

284 

15.6 

3 

39.7 

1362 

74.8 

i 
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I 

I 


The  percent  of  total  variance  is  also  provided  in  this 
table.  Clearly  the  high  frequency  terms  dominate  the 
spectrum. 

The  correlations  required  in  the  Section  VI  study 
can  be  calculated  with  Table  K-I  data  and  equations  (K-8) 
through  (K-15) .  Other  than  the  model  parameters  above 
ijj,  h,  and  h’  are  the  only  required  inputs.  Appendix  F 
explains  the  calculation  of  i|i  and  the  transformation  of 
the  Q-matrix  from  the  transverse-inplane-down  coordinates 
to  the  east-north-up  coordinates  of  the  navigation  error 
propagation  model. 
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APPENDIX  L 


Gravity  Disturbance  Statistical  Models  Comparison 

This  appendix  presents  a  graphical  view  of  the  auto- 
and  cross-correlations  of  gravity  disturbance  components. 

The  correlations  are  calculated  using  three  different  func¬ 
tional  forms « 

a.  Linear  State  Space.  This  shaping  filter  model  is 
presented  in  Appendix  F  with  additional  discussion 
in  Sections  IV  and  VI. 

b.  Anomaly  Degree  Variance.  This  Legendre  polynomial 
form  of  representation  is  presented  in  Appendix  J 
with  additional  discussion  in  Sections  V  and  VI. 

c.  Attenuated  White  Noise.  This  model,  based  on  a 
subterranean  white  noise  anomalous  potential,  is 
presented  in  Appendix  K  with  additional  comments 
in  Section  VI. 

The  altitude  level  was  zero  for  this  comparison.  The  coordi¬ 
nate  frame  in  which  gravity  disturbances  are  assumed  to  be 
expressed  is  a  local-level,  transverse-inplane-down 
(crossrange-downrange-down)  frame  used  in  References  19, 

21,  and  23.  Since  all  three  models  are  isotropic,  the 
transverse  component  cross-correlations  with  both  inplane 
and  down  components  are  zero,  hence  are  not  plotted. 


Figure  L-l .  Transverse  (Crossrange)  Gravity  Disturbance 

Autocorrelations 
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Figure  L-3-  Inplane  (Downrange)  and  Down  Gravity- 


Disturbances  Cross-Correlations 


260 


008T  0001  00>T  008t  0001  000  000  00*  008  0  008- 

(jPBSUl)  tlOH«iajJOQ 


Figure  L-4.  Gravity  Anomaly  (Down  Disturbance) 
Autocorrelations 
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